/*************************************************************************
 *
 *   
 *
 * 
 *************************************************************************/

#include "PrimaryGeneratorAction.hh"

using namespace std;
using namespace CLHEP;
G4int evctr = 0;
G4int transctr = 0;
G4int betactr = 0;


double betadist(double Q, double e3, double sc){  
  double m3 = 0.511;
  if ( e3 < 0 || e3>Q) 
    return 0;
  e3+=m3;
  double p3 = sqrt(e3*e3 - m3*m3);
  double fval=(Q+m3-e3)*(Q+m3-e3)*p3*e3/sc;
  return fval;
}


PrimaryGeneratorAction::PrimaryGeneratorAction(HistoManager* HM, Data* Data)
  :histoManager(HM), data(Data)
{
  G4int i,j,index;
  data->EventNumber = 0;
  data->DetectedEventNumber = 0;
  data->BeamParticles = 0;
  
  kinematics = NULL;
  ThetaCmMin = NULL;
  ThetaCmMax = NULL;
  MaxAngularDistribution = NULL;
  AngularDistribution = NULL;
  RutherfordFactor = NULL;
  E_cm = NULL;
  CrossSectionIntegral = NULL;
  ICrossSectionIntegral = NULL;
  MaxOpticalPotential = NULL;
  OpticalPotential = NULL;
  ArealDensity = NULL;
  ScatteringProbability = NULL;
  level = NULL;
  feeding = NULL;
  
  G4int n_particle = 1;
  
  particleGun1 = new G4ParticleGun(n_particle);
  particleGun2 = new G4ParticleGun(n_particle);
  particleGun3 = new G4ParticleGun(n_particle);
  particleGun4 = new G4ParticleGun(n_particle);
  

  if(!data->Efficiency && !data->AlphaSource && !data->QuadAlphaSource && !data->BetaEnergy)
    {
      //setup instance of Kinematic-Class
      kinematics = new Kinematic(data->Projectile, data->TargetMaterial, data->settings->TargetThickness/(mg/cm2));
      
      //cout << "projectile " << data->Projectile->Name() << " target " << data->TargetMaterial->Name() << endl;
      //cout << "target thick " << data->settings->TargetThickness/(mg/cm2) << " beam energy " << data->BeamEnergy/MeV << " step size " << data->ThicknessStepSize/(mg/cm2) << endl;

      //create necessary splines (Energy vs. X, Rutherford-Cross-Section vs. X)
      EVsX = kinematics->EnergyVsThickness(data->BeamEnergy/MeV,data->ThicknessStepSize/(mg/cm2));
      
      //for(int i=0;i<10;i++)
      //	cout << i << " thickness " << i/10. << " Energy " << EVsX->Eval(i/10.) << endl;
      
      RVsX = kinematics->RutherfordVsThickness(data->BeamEnergy/MeV,data->ThicknessStepSize/(mg/cm2),true);
      MaxRutherford = RVsX->Eval(data->settings->TargetThickness/(mg/cm2));
      
      if(!data->IsotropeParticleDistribution)
	{
	  //estimate minimum cm-angle for recoils to be detected
	  double energy[2],ThetaLab[4];
	  energy[0] = data->BeamEnergy;
	  energy[1] = EVsX->Eval(data->settings->TargetThickness/(mg/cm2));

	  //set the angles covered by the detectors in the lab-frame
	  SetThetaLab(ThetaLab);

	  //if(data->VerboseLevel > 3)
	  //  {
	      cout<<"Lab. Angels: "<<ThetaLab[0]/degree<<" - "<<ThetaLab[1]/degree<<", "<<ThetaLab[2]/degree<<" - "<<ThetaLab[3]/degree<<endl;
	      //   }

	  //calculate minimal theta_cm from lab-angles
//	  MinimalThetaCm(ThetaLab, energy);
//	  MaximalThetaCm(ThetaLab, energy);
	  if(!data->Coulex){
	    IntThetaCmMin =0;
	    IntThetaCmMax =180;

	  }
  		
     ThetaCmMin = new double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];
     ThetaCmMax = new double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];
	  
	  for (i=0;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++){
			ThetaCmMin[i] = data->ThetaCmMin;
			ThetaCmMax[i] = data->ThetaCmMax;
	  }

	  //verbose output of results
	  if(data->VerboseLevel > -1){//coutcheck wimmer
	    for(i=0;i<data->NumberOfLevel;i++){
	      cout<<"Using: Transfer reaction of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->Target->A()<<data->Target->Name()<<", E_ex = "<<data->LevelEnergy[i]/MeV<<" MeV: theta_cm_min = "<<ThetaCmMin[i]/degree<<" degree: theta_cm_max = "<<ThetaCmMax[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" degree)"<<endl;
	    }
	    if(!data->Coulex){
	      for(i=data->NumberOfLevel;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++){
		cout<<"Using: Elastic Scattering of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->A()<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->Name()<<": theta_cm_min = "<<ThetaCmMin[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	      }
	    }
	  }
	  //angular distributions
	  MaxAngularDistribution = new G4double[data->NumberOfLevel];
	  AngularDistribution = new TGraph*[data->NumberOfLevel];
	  RutherfordFactor = new G4double[data->TargetMaterial->NumberOfElements()];
	  E_cm = new G4double[data->TargetMaterial->NumberOfElements()];
	  if(!data->Coulex){
	    CrossSectionIntegral = new G4double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];
	    ICrossSectionIntegral = new G4double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];
	  }
	  else{
	    CrossSectionIntegral = new G4double[data->NumberOfLevel];
	    ICrossSectionIntegral = new G4double[data->NumberOfLevel];
	  }
	  
	  for(i=0;i<data->TargetMaterial->NumberOfElements();i++)
	    {
	      RutherfordFactor[i] = pow(data->TargetMaterial->GetElement(i)->Z()*data->Projectile->Z()*pow(eplus/coulomb,2)/(16*PI*epsilon0/(farad/meter)),2)*pow(meter*joule,2);

	      E_cm[i] = ((sqrt(2.*EVsX->Eval(0.)*data->TargetMaterial->GetElement(i)->Mass()+pow(data->Projectile->Mass()+data->TargetMaterial->GetElement(i)->Mass(),2)) +
			  sqrt(2.*EVsX->Eval(data->settings->TargetThickness/(mg/cm2))*data->TargetMaterial->GetElement(i)->Mass()+pow(data->Projectile->Mass()+data->TargetMaterial->GetElement(i)->Mass(),2)))/2.
			 -data->Projectile->Mass()-data->TargetMaterial->GetElement(i)->Mass())*MeV;
	      if(data->VerboseLevel > 1)
		{
		  cout<<data->TargetMaterial->GetElement(i)->A()<<data->TargetMaterial->GetElement(i)->Name()<<" RutherfordFactor = "<<RutherfordFactor[i]/(pow(MeV,2)*millibarn)<<" MeV^2*mb, E_cm = "<<E_cm[i]/MeV<<" MeV, E_lab = "<<EVsX->Eval(0.)/MeV<<" - "<<EVsX->Eval(data->settings->TargetThickness/(mg/cm2))/MeV<<" MeV"<<endl;
		  //	      cout<<"RutherfordFactor = pow("<<data->TargetMaterial->GetElement(i)->Z()<<"*"<<data->Projectile->Z()<<"*pow("<<eplus<<",2)/(16*"<<PI<<"*"<<epsilon0<<"),2)*1e25 = pow("<<data->TargetMaterial->GetElement(i)->Z()*data->Projectile->Z()<<"*"<<pow(eplus,2)<<"/"<<(16*PI*epsilon0)<<",2) = pow("<<data->TargetMaterial->GetElement(i)->Z()*data->Projectile->Z()*pow(eplus,2)/(16*PI*epsilon0)<<",2)"<<endl;
		  //	      cout<<"1/(16*PI*epsilon0)^2 = "<<1/pow(16*PI*epsilon0,2)/(pow(meter*joule,2)/pow(coulomb,4))<<" m^2*J^2/C^4, (Z1*Z2)^2 = "<<pow((double)(data->TargetMaterial->GetElement(i)->Z()*data->Projectile->Z()),2)<<", e^4 = "<<pow(eplus/coulomb,4)<<" C^4"<<endl;
		  //	      cout<<"eplus = "<<eplus/coulomb<<" C = "<<eplus<<endl;
		  //	      cout<<"epsilon0 = "<<epsilon0/(farad/meter)<<" farad/meter = "<<epsilon0<<endl;
		  //	      cout<<"RutherfordFactor = "<<RutherfordFactor[i]/(pow(MeV,2)*millibarn)<<" MeV^2*mb = "<<RutherfordFactor[i]<<endl;
		}
	    }
	  CreateSplines(data->AngularDistributionFile);
	  if(!data->Coulex){	    
	    MaxOpticalPotential = new G4double[data->TargetMaterial->NumberOfElements()];
	    OpticalPotential = new TGraph*[data->TargetMaterial->NumberOfElements()];
	    
	    if(data->OpticalPotential)
	      CreateOpticalPotential(data->OpticalPotentialFile);
	    else
	      CreateOpticalPotential("Rutherford");
	    
	    //check integrals
	    for(i=0;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++){
	      //if cross section integral smaller 0 => integration failed => exit
	      if(CrossSectionIntegral[i] < 0){
		cout<<"Warning, exiting because cross section integral "<<i<<" is negative ("<<CrossSectionIntegral[i]<<")"<<endl;
		
		exit(1);
	      }
	    }
	    
	  //increase transfer cross sections
	    if(data->TransferCrossSection < -1){
	      for(index=0;index<data->NumberOfLevel;index++){
		if(data->VerboseLevel > 3){
		  cout<<"multiplying CrossSectionIntegral["<<index<<"] = "<<CrossSectionIntegral[index]<<" with "<<-data->TransferCrossSection;
		}
		CrossSectionIntegral[index] *= -data->TransferCrossSection;
		if(data->VerboseLevel > 3){
		  cout<<" => "<<CrossSectionIntegral[index]<<endl;
		}
	      }
	    }
	  }
	  //calculate areal density
	  ArealDensity = new G4double[data->TargetMaterial->NumberOfElements()];
	  G4double TargetMass = data->TargetMaterial->Mass()/amu_c2*g;
	  if(data->VerboseLevel > 3){
	    cout<<"TargetMass = "<<TargetMass/mg<<" mg/mol ("<<data->TargetMaterial->Mass()<<", "<<amu_c2<<")"<<endl;
	  }
	  for(i=0;i<data->TargetMaterial->NumberOfElements();i++){
	    //is ->Fraction() right??
	    ArealDensity[i] = data->TargetMaterial->Fraction(i)*data->settings->TargetThickness/TargetMass*Avogadro;
	    if(data->VerboseLevel > 3){
	      cout<<"ArealDensity["<<i<<"] ("<<data->TargetMaterial->GetElement(i)->A()<<data->TargetMaterial->GetElement(i)->Name()<<") = "<<data->TargetMaterial->Fraction(i)<<"*"<<data->settings->TargetThickness/(mg/cm2)<<" mg/cm2/"<<TargetMass/mg<<" mg*"<<Avogadro<<" = "<<ArealDensity[i]/(1/cm2)<<"/cm2"<<endl;
	    }
	  }
	  
	  //sum up cross sections 
	  if(!data->Coulex){
	    ScatteringProbability = new G4double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];
	    ScatteringProbability[0] = CrossSectionIntegral[0]*ArealDensity[data->TargetMaterial->NumberOfElements()-1];
	    if(data->VerboseLevel > 3){      
	      cout<<"setting ScatteringProbability[0] to "<<ScatteringProbability[0]<<" (CrossSectionIntegral = "<<CrossSectionIntegral[0]/millibarn<<" mb, ArealDensity = "<<ArealDensity[data->TargetMaterial->NumberOfElements()-1]*millibarn<<" /mb)"<<endl;
	      
	    }
	    for(index=1;index<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();index++){
	      if(index < data->NumberOfLevel){
		ScatteringProbability[index] = ScatteringProbability[index-1] + ArealDensity[data->TargetMaterial->NumberOfElements()-1]*CrossSectionIntegral[index];
	      }
	      else{
		ScatteringProbability[index] = ScatteringProbability[index-1] + ArealDensity[index-data->NumberOfLevel]*CrossSectionIntegral[index];
	      }
	      if(data->VerboseLevel > 3){
		cout<<"setting ScatteringProbability["<<index<<"] to "<<ScatteringProbability[index]<<" (CrossSectionIntegral = "<<CrossSectionIntegral[index]/millibarn<<" mb";
		if(index >= data->NumberOfLevel){
		  cout<<", ArealDensity["<<index-data->NumberOfLevel<<"] = "<<ArealDensity[index-data->NumberOfLevel]*millibarn<<" /mb";
		}
		cout<<")"<<endl;
	      }
	    }
	  }
	  else{
	    //cout << "data->NumberOfLevel " << data->NumberOfLevel << " data->TargetMaterial->NumberOfElements() " << data->TargetMaterial->NumberOfElements() << endl;
	    ScatteringProbability = new G4double[data->NumberOfLevel+1];
	    //cout << "ArealDensity["<<data->TargetMaterial->NumberOfElements()-1<<"] = " << ArealDensity[data->TargetMaterial->NumberOfElements()-1] << endl;
	    //cout << "CrossSectionIntegral[0] " << CrossSectionIntegral[0] << endl;
	    ScatteringProbability[0] = ArealDensity[data->TargetMaterial->NumberOfElements()-1]*CrossSectionIntegral[0];
	    ScatteringProbability[data->NumberOfLevel] = ArealDensity[data->TargetMaterial->NumberOfElements()-1]*ICrossSectionIntegral[0];
	    //cout << "ScatteringProbability[0] " << ScatteringProbability[0] << endl;
	      if(data->VerboseLevel > 3){
		cout<<"setting ScatteringProbability[0] to "<<ScatteringProbability[0]<<" (CrossSectionIntegral = "<<CrossSectionIntegral[0]/millibarn<<" mb" <<" ICrossSectionIntegral = "<<ICrossSectionIntegral[0]/millibarn<<" mb)"<< endl;
	      }

	    for(index=1;index<data->NumberOfLevel;index++){
	      ScatteringProbability[index] = ScatteringProbability[index-1] + ArealDensity[data->TargetMaterial->NumberOfElements()-1]*CrossSectionIntegral[index];
	      ScatteringProbability[data->NumberOfLevel] += ArealDensity[data->TargetMaterial->NumberOfElements()-1]*ICrossSectionIntegral[index];
	      //cout << "ArealDensity["<<data->TargetMaterial->NumberOfElements()-1<<"] = " << ArealDensity[data->TargetMaterial->NumberOfElements()-1] << endl;
	      //cout << "CrossSectionIntegral["<<index<<"] " << CrossSectionIntegral[index] << endl;

	      if(data->VerboseLevel > 3){
		cout<<"setting ScatteringProbability["<<index<<"] to "<<ScatteringProbability[index]<<" (CrossSectionIntegral = "<<CrossSectionIntegral[index]/millibarn<<" mb" <<" ICrossSectionIntegral = "<<ICrossSectionIntegral[index]/millibarn<<" mb)" << endl;
	      }
	    }
	    cout<<"setting ScatteringProbability["<<data->NumberOfLevel<<"] = " << ScatteringProbability[data->NumberOfLevel] << endl;
	  }
	  
	}//if(!data->IsotropeParticleDistribution)
      //cout << "ScatteringProbability[data->NumberOfLevel-1] " << ScatteringProbability[data->NumberOfLevel-1] << endl;
      
      if(data->TestFlag)
	{
	  G4String filename = "test_VsX_" + G4UIcommand::ConvertToString(data->NumberOfThicknessSteps) + ".dat";
	  file.open(filename.data());
	  
	  for(G4double i=0;i<=data->settings->TargetThickness;i+=data->ThicknessStepSize)
	    {
	      file<<i/(mg/cm2)<<" "<<EVsX->Eval(i/(mg/cm2))<<" "<<RVsX->Eval(i/(mg/cm2))<<endl;
	    }
	  file.close();
	  
	  filename = data->ScratchDirectory+"/test_AngularDistribution.dat";
	  file.open(filename.data());
	}
      
      //create levelstructure and so on for gamma emission
      if(data->Gammas)
	{
	  level = new MiniBallSourceLevel*[data->NumberOfLevel];
	  feeding = new MiniBallSourceLevel();
	  for(G4int level_id=0; level_id < data->NumberOfLevel; data->LevelID[level_id] = level_id++)
	    {
	      level[level_id] = new MiniBallSourceLevel(data->LevelEnergy[level_id], data->LevelSpin[level_id], data->LevelParity[level_id]);
	      level_structure.push_back(level[level_id]);
	      feeding->decay.push_back(level_id);
	      feeding->prob.push_back(data->FeedingProbability[level_id]);
	      if(data->VerboseLevel > 4)//?
		{
		  cout<<"created level "<<level_id<<" ("<<data->LevelEnergy[level_id]<<", "<<data->LevelSpin[level_id]<<", "<<data->LevelParity[level_id]<<") with FeedingProbability "<<data->FeedingProbability[level_id]<<endl;
		}
	      for(G4int decay=0; decay < data->NumberOfDecays[level_id]; decay++)
		{
		  level[level_id]->AddDecay(data->LevelID[data->DecayLevel[level_id][decay]], data->DecayProbability[level_id][decay], data->DecayType[level_id][decay], data->DecayDelta[level_id][decay]);
		  if(data->VerboseLevel > 4)//?
		    {
		      cout<<"added decay "<<decay<<" ("<<data->LevelID[data->DecayLevel[level_id][decay]]<<", "<<data->DecayProbability[level_id][decay]<<", "<<data->DecayType[level_id][decay]<<", "<<data->DecayDelta[level_id][decay]<<")"<<endl;
		    }
		}
	    }
	  
	  // set internal conversion coefficients
	  //  SetInternalConversionCoeffE2(1.1716, 121.7825*keV);
	  
	  // correct relative decay probabilities by conversion coefficients
	  //InternalConversionCorrectProbabilities();
	  
	} //if(data->Gammas)
    }//if(!data->Efficiency)
  else
    {
      EfficiencyEnergy[0] = 121.7817;
      EfficiencyEnergy[1] = 244.6975;
      EfficiencyEnergy[2] = 344.2785;
      EfficiencyEnergy[3] = 778.9045;
      EfficiencyEnergy[4] = 964.1;
      EfficiencyEnergy[5] = 1085.836;
      EfficiencyEnergy[6] = 1112.074;
      EfficiencyEnergy[7] = 1408.011;
    }

  sum = 0;
  part_sum = new G4double[data->TargetMaterial->NumberOfElements()];
  for(i=0;i<data->TargetMaterial->NumberOfElements();i++)
    {
      sum += pow((double)data->TargetMaterial->GetElement(i)->Z(),2)/(pow(data->Projectile->A(),2)+pow(data->TargetMaterial->GetElement(i)->A(),2)+2.*data->TargetMaterial->GetElement(i)->A()*(data->Projectile->A()+data->BeamEnergy));
      part_sum[i] = 0;
      for(j=0;j<=i;j++)
	{
	  part_sum[i] += pow((double)data->TargetMaterial->GetElement(j)->Z(),2)/(pow(data->Projectile->A(),2)+pow(data->TargetMaterial->GetElement(j)->A(),2)+2.*data->TargetMaterial->GetElement(j)->A()*(data->Projectile->A()+data->BeamEnergy));
	}
    }

  particleTable = G4ParticleTable::GetParticleTable();
}

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  file.close();

  if(data->VerboseLevel > 1)
    {
      cout<<"starting PrimaryGeneratorAction Destructor"<<endl;
    }

  for(size_t q=0; q<level_structure.size(); q++)
    delete level_structure[q];

  if(data->VerboseLevel > 1)
    {
      cout<<"deleted level_structure"<<endl;
    }

  delete particleGun1;
  delete particleGun2;
  delete particleGun3;
  delete particleGun4;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted particleGuns"<<endl;
    }

  if(kinematics != NULL)
    delete kinematics;

  if(data->VerboseLevel > 1)
    {
      cout<<"deleted kinematics"<<endl;
    }

  if(ThetaCmMin != NULL)
    delete ThetaCmMin;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted ThetaCmMin"<<endl;
    }

  if(MaxAngularDistribution != NULL)
    delete[] MaxAngularDistribution;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted MaxAngularDistribution"<<endl;
    }

  if(AngularDistribution != NULL)
    delete[] AngularDistribution;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted AngularDistribution"<<endl;
    }

  if(RutherfordFactor != NULL)
    delete[] RutherfordFactor;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted RutherfordFactor"<<endl;
    }

  if(E_cm != NULL)
    delete[] E_cm;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted E_cm"<<endl;
    }

  if(CrossSectionIntegral != NULL)
    delete[] CrossSectionIntegral;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted CrossSectionIntegral"<<endl;
    }
  if(ICrossSectionIntegral != NULL)
    delete[] ICrossSectionIntegral;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted ICrossSectionIntegral"<<endl;
    }

  if(MaxOpticalPotential != NULL)
    delete[] MaxOpticalPotential;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted MaxOpticalPotential"<<endl;
    }

  if(OpticalPotential != NULL)
    delete[] OpticalPotential;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted OpticalPotential"<<endl;
    }

  if(ArealDensity != NULL)
    delete[] ArealDensity;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted ArealDensity"<<endl;
    }

  if(ScatteringProbability != NULL)
    delete[] ScatteringProbability;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted ScatteringProbability"<<endl;
    }


//  for(int i=0;i<data->NumberOfLevel;i++)
//    {
//      //cout<<"deleting level["<<i<<"] = "<<level[i]<<endl;
//      if(level[i] != NULL)
//	delete level[i];
//    }
//  if(data->VerboseLevel > 1)
//    {
//      cout<<"deleted levels"<<endl;
//    }

  if(level != NULL)
    delete[] level;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted level"<<endl;
    }

  if(feeding != NULL)
    delete feeding;
  if(data->VerboseLevel > 1)
    {
      cout<<"deleted feeding"<<endl;
    }

  delete[] part_sum;
  if(data->VerboseLevel > 1)
    {
      cout<<"DONE!"<<endl;
    }
}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
//  cout << "generate" << endl;
  
  //get point of reaction (from rutherford in z-direction and gaussian in x-y-direction)
  G4int gamma_nb;
  G4int i = 0;
  G4int reaction = 0;
  G4double reaction_x = 0.;
  G4double reaction_y = 0.;
  G4double reaction_z = 0.;
  G4double reaction_tmp = 0.;
  G4double theta_x = 0.;
  G4double theta_y = 0.;
  G4double theta_tmp = 0.;
  G4double tmp;
  G4double eex = 0;
  G4double phi;
  G4double phiN;
  G4double testv;
  //  G4bool transfer_reaction = true;

  G4double** energy;
  G4double** theta;

  //  G4double* blevelenergy; //kwi
  //  G4double* bpopulation;

  G4double reaction_energy = 0.;

  G4ThreeVector direction1;
  G4ThreeVector direction2;
  G4ThreeVector direction3;
  G4ThreeVector direction4;

  Element* rc = NULL;

  //kw add beta source
  if(data->AddBeta){
    testv = RandFlat::shoot(0.,1.);
    if(data->VerboseLevel > 4)
      cout << "test value = " << testv << endl; 
    if(testv < data->RatioBeta/(1+data->RatioBeta))
      data->BetaEnergy = data->AddBeta;
    else
      data->BetaEnergy = 0.;
  }
  // end add beta source

  //select random x and y position on a disk with diameter data->BeamWidth
  do {
    reaction_x = RandFlat::shoot(-data->BeamWidth/2.,data->BeamWidth/2.)*mm;
    reaction_y = RandFlat::shoot(-data->BeamWidth/2.,data->BeamWidth/2.)*mm;
    reaction_tmp = RandFlat::shoot(0.,data->BeamWidth/2.)*mm;
  } while(sqrt(pow(reaction_x,2)+pow(reaction_y,2)) > data->BeamWidth/2.);
//  } while(sqrt(pow(reaction_x,2)+pow(reaction_y,2)) > reaction_tmp);
/*  
  reaction_x = RandGauss::shoot(0,data->BeamWidth)*mm;
  reaction_y = RandGauss::shoot(0,data->BeamWidth)*mm;
*/
  data->ReactionX = reaction_x;
  data->ReactionY = reaction_y;

  //select random z-position
  if(!data->Efficiency && !data->DeltaE && !data->AlphaSource && !data->QuadAlphaSource && !data->BetaEnergy)
    {
/*      //quick and dirty method to generate a probability distribution in z-direction that matches the rutherford cross section
      do {
	reaction_z = RandFlat::shoot(0.,data->settings->TargetThickness/(mg/cm2))*mg/cm2;
	tmp = RandFlat::shoot(0.,MaxRutherford);
	if(data->VerboseLevel > 4)
	  {
	    i++;
	    cout<<"rutherford "<<i<<": "<<tmp<<"("<<0<<"-"<<MaxRutherford<<"), RVsX->Eval("<<reaction_z/(mg/cm2)<<") = "<<RVsX->Eval(reaction_z/(mg/cm2))<<" ('"<<0<<"-"<<data->settings->TargetThickness/(mg/cm2)<<"')"<<endl;
	  }
      } while(tmp > RVsX->Eval(reaction_z/(mg/cm2)));
      
      if(data->TestFlag)
	{
	  ofstream testfile;
	  G4String filename = data->ScratchDirectory+"/test_RVsX.dat";
	  testfile.open(filename.data(),ios::app);
	  
	  testfile<<tmp<<" "<<RVsX->Eval(reaction_z/(mg/cm2))<<" "<<reaction_z/(mg/cm2)<<endl;
	  
	  testfile.close();
	}
      
      data->TargetZ = reaction_z;
    }
  else if(!data->AlphaSource && !data->QuadAlphaSource) //beta source not on surface of target.
    {
  */    //flat distribution along depth of target
      //reaction_z = RandFlat::shoot(-data->settings->TargetThickness/(2.*data->settings->TargetMaterialDensity)/micrometer, data->settings->TargetThickness/(2.*data->settings->TargetMaterialDensity)/micrometer)*micrometer;
      reaction_z = RandFlat::shoot(0., data->settings->TargetThickness);
//      reaction_z = data->settings->TargetThickness/2.;
		data->ReactionZ = reaction_z;
	 
//	 	cout << "Reaction z: " << reaction_z/(mg/cm2) << endl;
	 }
  else
    {
      //alpha source: only from the surface of the 'target' (i.e. the source)
      tmp = RandFlat::shoot(0.,2.);
      if(tmp < 1.)
	{
	  reaction_x = -data->settings->TargetShift[0];
	  reaction_y = -data->settings->TargetShift[1];
	  reaction_z = -data->settings->TargetShift[2];
	  if(data->VerboseLevel > 1)
	    {
	      cout<<"tmp ("<<tmp<<") < 1. = "<<(tmp < 1.)<<" => reaction_z = "<<reaction_z/micrometer<<" um"<<endl;
	    }
	}
      else
	{
	  reaction_x = +data->settings->TargetShift[0];
	  reaction_y = +data->settings->TargetShift[1];
	  reaction_z = +data->settings->TargetShift[2];
	  if(data->VerboseLevel > 1)
	    {
	      cout<<"tmp ("<<tmp<<") < 1. = "<<(tmp < 1.)<<" => reaction_z = "<<reaction_z/micrometer<<" um"<<endl;
	    }
	  reaction_x *= -1;
	}
    }
  /*
  if(!data->AlphaSource  && !data->QuadAlphaSource){
    cout << " reaction_z = " <<reaction_z/micrometer<<" um " << endl;
    reaction_z+=data->settings->TargetShift;
    if(data->VerboseLevel > -1){
      cout << " reaction_z = " <<reaction_z/micrometer<<" um shifted by TargetShift = "<< data->settings->TargetShift/mm << " mm" <<endl;
    }
  }
  //doesn't work since reaction_z is in mg/cm2 ???
  */
  if(!data->Efficiency && !data->DeltaE && !data->AlphaSource  && !data->QuadAlphaSource && !data->BetaEnergy)
    {
      reaction_energy = EVsX->Eval(reaction_z/(mg/cm2))*MeV;
      data->ReactionEnergy = reaction_energy;
//      cout << "Beam energy " << data->BeamEnergy/MeV << endl;
//		cout << "reaction_z/(mg/cm2) " << reaction_z/(mg/cm2)  << " reaction_energy " << reaction_energy<< endl;
      if(data->VerboseLevel > 4){
	evctr++;
	transctr++;
	cout << "total: " << evctr << " transfer: " << transctr << endl; 
      }
      //choose reaction type
      reaction = ReactionType();
      //select excitation energy for transfer reactions
      eex = ExcitationEnergy(reaction);
      //generate theta_cm-distribution, minimum angle ThetaCmMin[] degree
      ThetaCm(reaction, eex);
      
      //select element from target material as elastic scattering partner
      if(data->Transfer == 0)
	ElasticTarget(reaction, rc);

      //set phi
      phi = RandFlat::shoot(-PI,PI)*radian;
      
      G4ParticleDefinition* ejectile = NULL;
      G4ParticleDefinition* recoil = NULL;
      G4ParticleDefinition* neutron = NULL;
      
      if(data->Transfer == 1) //transfer reaction
	{
	  //set particleGun1 to ejectile
	  if(data->VerboseLevel > 2)
	    cout<<"defining transfer ejectile ("<<flush<<data->Ejectile->Z()<<", "<<data->Ejectile->N()<<")"<<endl;

	  ejectile = Definition(data->Ejectile->Z(),data->Ejectile->N(),eex);
	  particleGun1->SetParticleDefinition(ejectile);
      
	  //set particleGun1 position and direction

	  //cout << "kinematics->orbits(data->Projectile,data->Target,data->Recoil,data->Ejectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,0,energy,theta)" <<endl;
	  //cout << "kinematics->orbits("<<data->Projectile<<","<<data->Target<<","<<data->Recoil<<","<<data->Ejectile<<","<<reaction_energy/MeV<<","<<data->Theta_cm/degree<<","<<data->Theta_cm/degree<<","<<1<<","<<eex/MeV<<","<<eex/MeV<<","<<1<<","<<0<<","<<energy<<","<<theta<<")" <<endl;

	  kinematics->orbits(data->Projectile,data->Target,data->Recoil,data->Ejectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,0,energy,theta);
      
	  //convert reaction_z from mg/cm2 to mm and shift so that target front at -thickness/2
	  reaction_z -= data->settings->TargetThickness/2.;
	  reaction_z /= data->settings->TargetMaterialDensity*10;


	  ///kwi_reaction_z
	  if(data->VerboseLevel > 1)
	    cout << " reaction_z = " <<reaction_z/micrometer<<" um " << endl;
	  reaction_x+=data->settings->TargetShift[0];
	  reaction_y+=data->settings->TargetShift[1];
	  reaction_z+=data->settings->TargetShift[2];
	  if(data->VerboseLevel > 1)
	    cout << " reaction_z = " <<reaction_z/micrometer<<" um shifted by TargetShift = "<< data->settings->TargetShift[2]/mm << " mm" <<endl;
	  

	  data->ReactionZ = reaction_z;
	  if(data->VerboseLevel > 2)
	    {
	      cout<<"reaction_z = "<<reaction_z/micrometer<<" um"<<endl;
	    }
      
	  particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
//	  cout << theta[0][0]*180./PI << "\t" << energy[0][0] << "\t" << endl;
	  energy[0][0] *= MeV;


    	 theta_x = RandFlat::shoot(0.,data->AngleWidth);
    	 theta_y = RandFlat::shoot(-180.,180.);
	 theta[0][0] += theta_x*cos(theta_y) ;
	 theta[0][0] += data->settings->BeamAngle*cos(phi - data->settings->BeamDirection);
	 theta[0][0] *= radian;
      
	  direction1.set(0,0,1);
	  direction1.setTheta(theta[0][0]);
	  direction1.setPhi(phi);
	  direction1.setMag(kinematics->Momentum(energy[0][0],data->Ejectile->Mass()));
	  particleGun1->SetParticleMomentumDirection(direction1);
	  //particleGun1->SetParticleMomentum(kinematics->Momentum(energy[0][0],data->Ejectile->Mass()));
	  if(data->OnePartOnly%10==2)
	    particleGun1->SetParticleEnergy(0.0);
	  else if(data->OnePartOnly==11)
	    particleGun1->SetParticleEnergy(100*MeV);
	  else 
	    particleGun1->SetParticleEnergy(energy[0][0]);
	  //cout << "momentum p " << kinematics->Momentum(energy[0][0],data->Ejectile->Mass()) << "\t p " << particleGun1->GetParticleMomentum()<< "\t e " << particleGun1->GetParticleEnergy() << endl; // momentum of ejectile
      
	  data->EjectileEnergy = energy[0][0];
	  data->EjectileTheta = theta[0][0];
	  data->EjectilePhi = phi;
	 
//	 cout << "E1: " << data->EjectileEnergy << endl;
      
	  delete[] energy[0];
	  delete[] energy;
	  delete[] theta[0];
	  delete[] theta;
      
	  //set particleGun2 to recoil
	  if(data->VerboseLevel > 2)
	    cout<<"defining transfer recoil ("<<flush<<data->Recoil->Z()<<", "<<data->Recoil->N()<<")"<<endl;

	  recoil = Definition(data->Recoil->Z(),data->Recoil->N(),0.0);

	  particleGun2->SetParticleDefinition(recoil);
      
	  //set particleGun2 position and direction
	  kinematics->orbits(data->Projectile,data->Target,data->Recoil,data->Ejectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,1,energy,theta);
      
	  particleGun2->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
	  //cout << theta[0][0]*180./PI << "\t" << energy[0][0] << endl;
	  energy[0][0] *= MeV;
  	  
    	 theta_x = RandFlat::shoot(0.,data->AngleWidth);
    	 theta_y = RandFlat::shoot(-180.,180.);
	 theta[0][0] += theta_x*cos(theta_y) ;
	  theta[0][0] += data->settings->BeamAngle*cos(phi+180.*degree - data->settings->BeamDirection);
	  theta[0][0] *= radian;
      
	  direction2.set(0,0,1);
	  direction2.setTheta(theta[0][0]);
	  direction2.setPhi(phi+180.*degree);
	  particleGun2->SetParticleMomentumDirection(direction2);
	  if(data->OnePartOnly%10==1)
	    particleGun2->SetParticleEnergy(0.0);
	  else if(data->OnePartOnly==12)
	    particleGun2->SetParticleEnergy(100*MeV);
	  else
	    particleGun2->SetParticleEnergy(energy[0][0]);
	  
	  data->RecoilEnergy = energy[0][0];
	  data->RecoilTheta = theta[0][0];
	  data->RecoilPhi = phi+180.*degree;
	 
//	 cout << "E2: " << data->RecoilEnergy << endl;

	  delete[] energy[0];
	  delete[] energy;
	  delete[] theta[0];
	  delete[] theta;
	  
/*	  //set particleGun4 to neutron
	  neutron = Definition(0,1,eex);
	  particleGun4->SetParticleDefinition(neutron);
      
	  kinematics->orbits(data->Projectile,data->Target,data->Recoil,data->Ejectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,1,energy,theta);
      
	  particleGun4->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
//	  cout << theta[0][0]*180./PI << "\t" << energy[0][0] << "\t" << endl;
	  energy[0][0] *= MeV;
	  theta[0][0] += data->settings->BeamAngle*cos(phiN - data->settings->BeamDirection);
	  theta[0][0] *= radian;
      
	  direction4.set(0,0,1);
	  direction4.setTheta(theta[0][0]);
	  direction4.setPhi(phiN);
	  direction4.setMag(kinematics->Momentum(energy[0][0],1.));
	  particleGun4->SetParticleMomentumDirection(direction4);
	  //particleGun1->SetParticleMomentum(kinematics->Momentum(energy[0][0],data->Ejectile->Mass()));
	  if(data->OnePartOnly%10==2)
	    particleGun4->SetParticleEnergy(0.0);
	  else if(data->OnePartOnly==11)
	    particleGun4->SetParticleEnergy(100*MeV);
	  else 
	    particleGun4->SetParticleEnergy(energy[0][0]);
	  //cout << "momentum p " << kinematics->Momentum(energy[0][0],data->Ejectile->Mass()) << "\t p " << particleGun1->GetParticleMomentum()<< "\t e " << particleGun1->GetParticleEnergy() << endl; // momentum of ejectile
      
	  data->EjectileEnergy = energy[0][0];
	  data->EjectileTheta = theta[0][0];
	  data->EjectilePhi = phi;
      
	  delete[] energy[0];
	  delete[] energy;
	  delete[] theta[0];
	  delete[] theta;
*/	}//if(data->Transfer == 1)
      else //elastic scattering => ejectile = projectile, recoil = (elastic)target
	{
	  //set particleGun1 to ejectile
	  if(data->VerboseLevel > 4)
	    cout<<"defining elastic ejectile ("<<flush<<data->Projectile->Z()<<", "<<data->Projectile->N()<<")"<<endl;

	  ejectile = Definition(data->Projectile->Z(),data->Projectile->N(),eex);

	  particleGun1->SetParticleDefinition(ejectile);
      
	  //set particleGun1 position and direction
	  kinematics->orbits(data->Projectile,rc,rc,data->Projectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,0,energy,theta);
      
	  //convert reaction_z from mg/cm2 to mm and shift so that target front at -thickness/2
	  reaction_z -= data->settings->TargetThickness/2.;
	  reaction_z /= data->settings->TargetMaterialDensity;

	  ///kwi_reaction_z
	  if(data->VerboseLevel > 1)
	    cout << " reaction_z = " <<reaction_z/micrometer<<" um " << endl;
	  reaction_x+=data->settings->TargetShift[0];
	  reaction_y+=data->settings->TargetShift[1];
	  reaction_z+=data->settings->TargetShift[2];
	  if(data->VerboseLevel > 1)
	    cout << " reaction_z = " <<reaction_z/micrometer<<" um shifted by TargetShift = "<< data->settings->TargetShift[2]/mm << " mm" <<endl;

	  data->ReactionZ = reaction_z;
	  if(data->VerboseLevel > 4)
	    {
	      cout<<"reaction_z = "<<reaction_z/micrometer<<" um"<<endl;
	    }
      
	  particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
	  energy[0][0] *= MeV;
  	  
    	 theta_x = RandFlat::shoot(0.,data->AngleWidth);
    	 theta_y = RandFlat::shoot(-180.,180.);
	 theta[0][0] += theta_x*cos(theta_y) ;
	  theta[0][0] += data->settings->BeamAngle*cos(phi - data->settings->BeamDirection);
	  theta[0][0] *= radian;
	  // testing output
	  if(data->VerboseLevel > 0){
	    cout << "data->Reaction " << data->Reaction << endl;
	    cout << "ejectile z " << rc->Z() << " n " << rc->N()<<endl;
	    cout << "theta cm " << data->Theta_cm << "lab " << theta[0][0] << endl;
	  }
	  direction1.set(1,1,1);
	  direction1.setTheta(theta[0][0]);
	  direction1.setPhi(phi);
	  particleGun1->SetParticleMomentumDirection(direction1);
      
	  particleGun1->SetParticleEnergy(energy[0][0]);
      
	  data->EjectileEnergy = energy[0][0];
	  data->EjectileTheta = theta[0][0];
	  data->EjectilePhi = phi;
      
	  delete[] energy[0];
	  delete[] energy;
	  delete[] theta[0];
	  delete[] theta;
      
	  //set particleGun2 to recoil

	  if(data->VerboseLevel > 4)
	    cout<<"defining elastic recoil "<<i<<" ("<<flush<<rc->Z()<<", "<<rc->N()<<")"<<endl;

	  recoil = Definition(rc->Z(),rc->N(),0.0);

	  particleGun2->SetParticleDefinition(recoil);
      
	  //set particleGun2 position and direction
	  kinematics->orbits(data->Projectile,rc,rc,data->Projectile,reaction_energy/MeV,data->Theta_cm/degree,data->Theta_cm/degree,1,eex/MeV,eex/MeV,1,1,energy,theta);
      
	  particleGun2->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
	  energy[0][0] *= MeV;
  	  
    	 theta_x = RandFlat::shoot(0.,data->AngleWidth);
    	 theta_y = RandFlat::shoot(-180.,180.);
	 theta[0][0] += theta_x*cos(theta_y) ;
	  theta[0][0] += data->settings->BeamAngle*cos(phi+180.*degree - data->settings->BeamDirection);
	  theta[0][0] *= radian;
      
	  direction2.set(0,0,1);
	  direction2.setTheta(theta[0][0]);
	  direction2.setPhi(phi+180.*degree);
	  particleGun2->SetParticleMomentumDirection(direction2);
      
	  particleGun2->SetParticleEnergy(energy[0][0]);
      
	  data->RecoilEnergy = energy[0][0];
	  data->RecoilTheta = theta[0][0];
	  data->RecoilPhi = phi+180.*degree;
      
	  delete[] energy[0];
	  delete[] energy;
	  delete[] theta[0];
	  delete[] theta;
	}//end of elastic scattering

      particleGun1->GeneratePrimaryVertex(anEvent);
      particleGun2->GeneratePrimaryVertex(anEvent); 
  //    particleGun4->GeneratePrimaryVertex(anEvent);     

      if(data->VerboseLevel > 4)
	{
	  cout<<"created event "<<data->EventNumber+1<<" with E = "<<data->ReactionEnergy/MeV<<" MeV, data->Theta_cm = "<<data->Theta_cm/degree<<" degree and Eex = "<<data->EexEnergy/MeV<<" MeV:"<<endl;
	    if(data->Transfer == 1)
	      {
		cout<<"Ejectile ("<<data->Ejectile->A()<<data->Ejectile->Name()<<"/"<<ejectile->GetParticleName()<<" from A: "<<data->Ejectile->Z()+data->Ejectile->N()<<", Z: "<<data->Ejectile->Z()<<"): \tE = "<<data->EjectileEnergy/MeV<<" MeV, Theta = "<<data->EjectileTheta/degree<<" degree and Phi = "<<data->EjectilePhi/degree<<" degree ("<<direction1.getTheta()/degree<<", "<<direction1.getPhi()/degree<<")"<<endl
		    <<"Recoil ("<<data->Recoil->A()<<data->Recoil->Name()<<"/"<<recoil->GetParticleName()<<" from A: "<<data->Recoil->Z()+data->Recoil->N()<<", Z: "<<data->Recoil->Z()<<"): \tE = "<<data->RecoilEnergy/MeV<<" MeV, Theta = "<<data->RecoilTheta/degree<<" degree and Phi = "<<data->RecoilPhi/degree<<" degree ("<<direction2.getTheta()/degree<<", "<<direction2.getPhi()/degree<<")"<<endl;
	      }
	    else
	      {
		cout<<"Ejectile ("<<data->Projectile->A()<<data->Projectile->Name()<<"/"<<ejectile->GetParticleName()<<" from A: "<<data->Projectile->Z()+data->Projectile->N()<<", Z: "<<data->Projectile->Z()<<"): \tE = "<<data->EjectileEnergy/MeV<<" MeV, Theta = "<<data->EjectileTheta/degree<<" degree and Phi = "<<data->EjectilePhi/degree<<" degree ("<<direction1.getTheta()/degree<<", "<<direction1.getPhi()/degree<<")"<<endl
		    <<"Target ("<<rc->A()<<rc->Name()<<"/"<<recoil->GetParticleName()<<" from A: "<<rc->Z()+rc->N()<<", Z: "<<rc->Z()<<"): \tE = "<<data->RecoilEnergy/MeV<<" MeV, Theta = "<<data->RecoilTheta/degree<<" degree and Phi = "<<data->RecoilPhi/degree<<" degree ("<<direction2.getTheta()/degree<<", "<<direction2.getPhi()/degree<<")"<<endl;
	      }
	}


      if(data->Gammas)
	{
	  G4LorentzVector particle, boosted_gamma;
	  if(data->Transfer == 1)
	    {
	      if(data->NumberOfDecays[reaction] > 0)
		{
		  if(data->VerboseLevel > 5)
		    {
		      cout<<"trying to create gammas:"<<endl;
		      for(gamma_nb=0;gamma_nb<data->NumberOfDecays[reaction];gamma_nb++)
			{
			  cout<<gamma_nb<<".: Level = "<<data->DecayLevel[reaction][gamma_nb]<<", Probability = "<<data->DecayProbability[reaction][gamma_nb]<<", Type = "<<data->DecayType[reaction][gamma_nb]<<", Delta = "<<data->DecayDelta[reaction][gamma_nb]<<endl;

			}
		    }

		 // GenerateDecay(reaction);
		  GenerateDecay();
		  
		  double p = 0.0569* 0.299;
		  if(data->Lifetime<0.1)
		    reaction_z = 0.;
		  else
		    reaction_z = RandExponential::shoot(p*data->Lifetime)*m;
		  //cout << " gamma from " << reaction_z << endl;
		  particleGun3->SetParticleDefinition(G4Gamma::GammaDefinition());
		  particleGun3->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
		  
		  if(data->VerboseLevel > 4)
		    {
		      cout<<"creating gammas:"<<endl;
		    }
		  for(gamma_nb=0; gamma_nb<GetNbOfGammas(); gamma_nb++)
		    {
		      if(data->IsoGamma){
			direction3.set(1,1,1);
			direction3.setPhi(RandFlat::shoot(-PI, PI)*radian);
			direction3.setTheta(RandFlat::shoot(0., PI)*radian);
			boosted_gamma.set(GetGammaEnergy(gamma_nb),direction3);
		      }
		      else
			boosted_gamma.set(GetGammaEnergy(gamma_nb),GetGammaDirection(gamma_nb));

		      //cout << boosted_gamma.phi() <<"\t"<< boosted_gamma.theta() << endl;
		      particle.set(data->EjectileEnergy+data->Ejectile->Mass(),direction1);
		      boosted_gamma.boost(particle.vect(),particle.beta());
		      particleGun3->SetParticleEnergy(boosted_gamma.e());

		      particleGun3->SetParticleMomentumDirection(boosted_gamma.vect());		      
		      /*
		      boosted_gamma.set(GetGammaEnergy(gamma_nb),GetGammaDirection(gamma_nb));
		      //cout << "energy " << data->EjectileEnergy+data->Ejectile->Mass() << "\t" << direction1.mag() << endl;
		      particle.set(data->EjectileEnergy+data->Ejectile->Mass(),direction1);
		      boosted_gamma.boost(particle.vect(),particle.beta());
		      particleGun3->SetParticleEnergy(boosted_gamma.e());
		      
		      if(data->IsoGamma){
			direction3.set(1,1,1);
			direction3.setPhi(RandFlat::shoot(-PI, PI)*radian);
			direction3.setTheta(RandFlat::shoot(0., PI)*radian);
			particleGun3->SetParticleMomentumDirection(direction3);
		      }
		      else
		        particleGun3->SetParticleMomentumDirection(boosted_gamma.vect());
		      
		      */
		      particleGun3->GeneratePrimaryVertex(anEvent);
		      //cout << "energy " << data->EjectileEnergy << "\ttheta " << direction1.getTheta()/degree << "\tbeta " << particle.beta() << "\tdopllershifted " << boosted_gamma.e() << "\tmag " << direction1.mag() << endl;
		      if(data->VerboseLevel > 4)
			{
			  cout<<"\t"<<gamma_nb<<".: E = "<<boosted_gamma.e()<<" (boosted from "<<GetGammaEnergy(gamma_nb)<<")"<<endl;
			}
		    }
		}
	    }
	  else
	    {
	    }
	}
    }//if(!data->Efficiency && !data->DeltaE)
  else if(data->DeltaE)
    {
      particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      
      if(data->VerboseLevel > 4)
	cout<<"shooting from ("<<reaction_x<<", "<<reaction_y<<", "<<reaction_z<<"): "<<flush;

      //same probability for protons, deuterons and tritons
      do {
	tmp = RandFlat::shoot(0.,3.);
      } while (tmp == 3);

      if(data->VerboseLevel > 4)
	cout<<tmp<<" = "<<flush;

      if(tmp < 1)
	{
	  particleGun1->SetParticleDefinition(G4Proton::ProtonDefinition());
	  if(data->VerboseLevel > 4)
	    cout<<"proton"<<endl;
	}
      else if(tmp < 2)
	{
	  particleGun1->SetParticleDefinition(G4Deuteron::DeuteronDefinition());
	  if(data->VerboseLevel > 4)
	    cout<<"deuteron"<<endl;
	}
      else
	{
	  particleGun1->SetParticleDefinition(G4Triton::TritonDefinition());
	  if(data->VerboseLevel > 4)
	    cout<<"triton"<<endl;
	}

      reaction = (int) tmp;

      //isotropic distribution
      data->Theta_cm = RandFlat::shoot(-1., 1.);
      data->Theta_cm = acos(data->Theta_cm)*radian;

      phi = RandFlat::shoot(-M_PI,M_PI)*radian;

      if(data->VerboseLevel > 4)
	cout<<"in direction "<<flush;

      direction1.set(0,0,1);
      direction1.setTheta(data->Theta_cm);
      direction1.setPhi(phi);
      particleGun1->SetParticleMomentumDirection(direction1);
      
      if(data->VerboseLevel > 4)
	cout<<"("<<data->Theta_cm<<", "<<phi<<")"<<endl;

      //flat energy dist.
      tmp = RandFlat::shoot(0.,30.)*MeV;
      particleGun1->SetParticleEnergy(tmp);

      data->RecoilTheta = data->Theta_cm;
      data->RecoilPhi = phi;
      data->RecoilEnergy = tmp;

      particleGun1->GeneratePrimaryVertex(anEvent);
    }//else if(data->DeltaE)
  else if(data->Efficiency&&data->Lifetime<0)
    {
      particleGun1->SetParticleDefinition(G4Gamma::GammaDefinition());
      particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));

      //calculate gamma energy
      do{i = ((G4int)RandFlat::shoot(0., 8.));} while(i>7);
      particleGun1->SetParticleEnergy(EfficiencyEnergy[i]*keV);

      data->ReactionEnergy = EfficiencyEnergy[i]*keV;

      //isotropic distribution
      data->Theta_cm = RandFlat::shoot(-1., 1.);
      data->Theta_cm = acos(data->Theta_cm)*radian;
      
      phi = RandFlat::shoot(-M_PI, M_PI)*radian; 
      
      direction1.set(0., 0., 1);
      direction1.setTheta(data->Theta_cm);
      direction1.setPhi(phi);
      
      particleGun1->SetParticleMomentumDirection(direction1);

      particleGun1->GeneratePrimaryVertex(anEvent);

      data->GammaCounter[i]++;
    }
  else if(data->Efficiency&&data->Lifetime>-0.5)
    {
      //c = 299792458m/s = 0.299 m/ns
      //beta = 0.0569 c
      // x = v*t = 0.0569* 0.299 * L meter
      double p = 0.0569* 0.299;
      if(data->Lifetime==0)
	reaction_z = 0.;
      else
	reaction_z = RandExponential::shoot(p*data->Lifetime)*m;

      //cout << "shooting from z = " << reaction_z << endl;
      particleGun1->SetParticleDefinition(G4Gamma::GammaDefinition());
      particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));

      particleGun1->SetParticleEnergy(886*keV);

      data->ReactionEnergy = 886*keV;

      //isotropic distribution
      data->Theta_cm = RandFlat::shoot(-1., 1.);
      data->Theta_cm = acos(data->Theta_cm)*radian;
      
      phi = RandFlat::shoot(-M_PI, M_PI)*radian; 
      
      direction1.set(0., 0., 1);
      direction1.setTheta(data->Theta_cm);
      direction1.setPhi(phi);
      
      particleGun1->SetParticleMomentumDirection(direction1);

      particleGun1->GeneratePrimaryVertex(anEvent);

      data->GammaCounter[i]++;
    }
  else if(data->AlphaSource)
    {
      //choose the emitted alpha      
      tmp = RandFlat::shoot(0.,3.);

      if(tmp < 1)//239Pu
	{
	  tmp = RandFlat::shoot(0., 99.82);//70.77+17.11+11.94

	  if(tmp < 70.77)
	    {
	      reaction_energy = 5156.59*keV;
	    }
	  else if(tmp < 70.77+17.11)
	    {
	      reaction_energy = 5144.3*keV;
	    }
	  else
	    {
	      reaction_energy = 5105.5*keV;
	    }
	}
      else if(tmp < 2)//241Am
	{
	  tmp = RandFlat::shoot(0., 97.9);//84.8+13.1

	  if(tmp < 84.8)
	    {
	      reaction_energy = 5485.56*keV;
	    }
	  else
	    {
	      reaction_energy = 5442.80*keV;
	    }
	}
      else//244Cm
	{
	  tmp = RandFlat::shoot(0., 100.);//76.4+23.6

	  if(tmp < 76.4)
	    {
	      reaction_energy = 5804.77*keV;
	    }
	  else
	    {
	      reaction_energy = 5762.64*keV;
	    }
	}

      particleGun1->SetParticleDefinition(Definition(2,2,0.0));
      particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      particleGun1->SetParticleEnergy(reaction_energy);

      //isotropic distribution, if reaction_z < 0 shoot only to negative z's (theta < 90) otherwise only to positive z's (theta > 90)
      if(reaction_z < 0)
	{
	  data->Theta_cm = RandFlat::shoot(-1., 0.);
	}
      else
	{
	  data->Theta_cm = RandFlat::shoot(0., 1.);
	}
      data->Theta_cm = acos(data->Theta_cm)*radian;

      phi = RandFlat::shoot(-M_PI,M_PI)*radian;

      if(data->VerboseLevel > 4)
	cout<<"in direction "<<flush;

      direction1.set(0,0,1);
      direction1.setTheta(data->Theta_cm);
      direction1.setPhi(phi);
      particleGun1->SetParticleMomentumDirection(direction1);

      particleGun1->GeneratePrimaryVertex(anEvent);

      data->ReactionZ = reaction_z;
      data->ReactionEnergy = reaction_energy;

      data->EjectileEnergy = tmp;
      data->EjectileTheta = data->Theta_cm;
      data->EjectilePhi = phi;
    }

  //kw quadruple alpha source 148Gd

  else if(data->QuadAlphaSource){
    //choose the emitted alpha      
    tmp = RandFlat::shoot(0.,4.);
    
    if(tmp < 1){ //239Pu
      tmp = RandFlat::shoot(0., 99.82);//70.77+17.11+11.94
      if(tmp < 70.77){
	reaction_energy = 5156.59*keV;
      }
      else if(tmp < 70.77+17.11){
	reaction_energy = 5144.3*keV;
      }
      else{
	reaction_energy = 5105.5*keV;
      }
    }
    else if(tmp < 2){ //241Am
      tmp = RandFlat::shoot(0., 97.9);//84.8+13.1
      if(tmp < 84.8){
	reaction_energy = 5485.56*keV;
      }
      else{
	reaction_energy = 5442.80*keV;
      }
    }
    else if(tmp < 3){ //244Cm
      tmp = RandFlat::shoot(0., 100.);//76.4+23.6
      if(tmp < 76.4){
	reaction_energy = 5804.77*keV;
      }
      else{
	reaction_energy = 5762.64*keV;
      }
    }
    else { //148Gd
      reaction_energy = 3182.69*keV;
    }
    
    particleGun1->SetParticleDefinition(Definition(2,2,0.0));
    particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
    particleGun1->SetParticleEnergy(reaction_energy);
    
      //isotropic distribution, if reaction_z < 0 shoot only to negative z's (theta < 90) otherwise only to positive z's (theta > 90)
    if(data->settings->TargetShift[2]<0.01){
      data->Theta_cm = RandFlat::shoot(-1., 1.);
    }
    else{
      if(reaction_z < 0){
	data->Theta_cm = RandFlat::shoot(-1., 0.);
      }
      else{
	data->Theta_cm = RandFlat::shoot(0., 1.);
      }
    }
    data->Theta_cm = acos(data->Theta_cm)*radian;
    
    phi = RandFlat::shoot(-M_PI,M_PI)*radian;
      
    if(data->VerboseLevel > 4)
      cout<<"in direction "<<flush;

    direction1.set(0,0,1);
    direction1.setTheta(data->Theta_cm);
    direction1.setPhi(phi);
    particleGun1->SetParticleMomentumDirection(direction1);
    
    particleGun1->GeneratePrimaryVertex(anEvent);
    
    data->ReactionZ = reaction_z;
    data->ReactionEnergy = reaction_energy;
    
    data->EjectileEnergy = tmp;
    data->EjectileTheta = data->Theta_cm;
    data->EjectilePhi = phi;
  }
  //kw beta source
  //  else if(data->BetaEnergy || (data->AddBeta && !transfer_reaction))
  else if(data->BetaEnergy) 
{

//    monoenergetic beta source 
//      reaction_energy = data->BetaEnergy*1000*keV;

//    flat distribution 0 - max MeV
//      reaction_energy = RandFlat::shoot(0., data->BetaEnergy)*1000*keV;

//    gaussian distribution mean = BetaEnergy/2, 
//      reaction_energy = RandGauss::shoot(data->BetaEnergy/2., data->BetaEnergy/3./ 2.35)*1000*keV;

//    Fit gauss to beta spectrum with fermi funtion, max 7 MeV, mean 3.235 MeV, sigma 1.5849 MeV
//    /home/kwimmer/transfer/simulation/beta/gausfit.ps
//      reaction_energy = RandGauss::shoot(3.235, 1.585)*1000*keV;


/*
//
//beta decay to excited levels

  int bnroflevels=3;
  blevelenergy = new double[bnroflevels];
  bpopulation = new double[bnroflevels];

  //mg30

  blevelenergy[0] = 0.0;
  blevelenergy[1] = 0.244;
  blevelenergy[2] = 0.688;
  
  bpopulation[0] = 0.3;
  bpopulation[1] = 0.7;
  bpopulation[2] = 0.1;
  
  double pop_tot=0;
  double qvalue=0*MeV;
  for(int i =0; i<bnroflevels; i++)
    pop_tot+=bpopulation[i];
  for(int i =0; i<bnroflevels; i++)
    bpopulation[i]/=pop_tot;

  tmp = RandFlat::shoot(0.,1.);

  for(int i =0; i<bnroflevels; i++){
    bpopulation[i]+=bpopulation[i-1];
    if(tmp<bpopulation[i]){      
      qvalue = data->BetaEnergy - blevelenergy[i];
      cout << i << " pop " << bpopulation[i] << " Q-value " << qvalue << endl;
      break;
    }
  }
  cout << endl;

//end beta decay to excited levels
*/
  double qvalue=0*MeV;
  tmp = RandFlat::shoot(0.,1.);
  for(int i =0; i<data->BetaNrOfLevels; i++){
    if(tmp < data->BetaPopulation[i]){      
      qvalue = data->BetaEnergy - data->BetaLevelEnergy[i];
      //cout << i << " pop " << data->BetaPopulation[i] << " levelen " << data->BetaLevelEnergy[i] << " Q-value " << qvalue << endl;
      break;
    }
  }
  double maxen=0;
  if(qvalue<1)
    maxen = 5.073333e-01 + 2.307576e-01*qvalue + 1.136364e-01*qvalue*qvalue;//low q values not linear
  else
    maxen = 3.285725e-01+4.923773e-01*qvalue;//fit /home/kwimmer/transfer/simulation/beta/maximum.ps
      
  double scaling = betadist(qvalue, maxen, 1); 

      do {
	tmp = RandFlat::shoot(0.,1.);
	reaction_energy = RandFlat::shoot(0.,qvalue);
      } while(tmp > betadist(qvalue, reaction_energy, scaling));
      

      if(data->VerboseLevel > 4){
	evctr++;
	betactr++;
	cout << "total: " << evctr << " beta: " << betactr << endl; 
      }

      //monoenergetic proton source for test 
      //reaction_energy = data->BetaEnergy*1000*keV;
      //particleGun1->SetParticleDefinition(G4Proton::ProtonDefinition ());

      particleGun1->SetParticleDefinition(G4Electron::ElectronDefinition ());
      particleGun1->SetParticlePosition(G4ThreeVector(reaction_x, reaction_y, reaction_z));
      particleGun1->SetParticleEnergy(reaction_energy);

      //isotropic distribution, if reaction_z < 0 shoot only to negative z's (theta < 90) otherwise only to positive z's (theta > 90)
      if(reaction_z < 0)
	{
	  data->Theta_cm = RandFlat::shoot(-1., 0.);
	}
      else
	{
	  data->Theta_cm = RandFlat::shoot(0., 1.);
	}
      data->Theta_cm = acos(data->Theta_cm)*radian;

      phi = RandFlat::shoot(-M_PI,M_PI)*radian;

      if(data->VerboseLevel > 4)
	cout<<"in direction "<<flush;

      direction1.set(0,0,1);
      direction1.setTheta(data->Theta_cm);
      direction1.setPhi(phi);
      particleGun1->SetParticleMomentumDirection(direction1);

      particleGun1->GeneratePrimaryVertex(anEvent);

      data->ReactionZ = reaction_z;
      data->ReactionEnergy = reaction_energy;

      data->EjectileEnergy = tmp;
      data->EjectileTheta = data->Theta_cm;
      data->EjectilePhi = phi;

    }
  ////////
  else
    {
      cerr<<"Error in "<<__PRETTY_FUNCTION__<<endl;
      exit(2);
    }

  data->Reaction = reaction;

  data->EventNumber++;

//	cout << "NEvent: "<< data->EventNumber << endl;

  G4int tenthousandth = (G4int)(data->EventNumber/10000);
  G4int hundredthousandth = (G4int)(data->EventNumber/100000);
  if( (data->EventNumber - hundredthousandth*100000)==0 )
    G4cout << data->EventNumber << G4endl;
  else if( (data->EventNumber - tenthousandth*10000)==0 )
    G4cout << "." << flush;
}

void PrimaryGeneratorAction::CreateSplines(G4String filename)
{
  G4int i,j,maxentries;
  G4int* entries;
  ifstream file;
  G4double** x;
  G4double** y;
  G4int first, last;
  G4int ifirst, ilast;
  double xtmp, ytmp;

  G4String stmp = data->SimulationDirectory+"/AngularDistribution.root";

  if(data->VerboseLevel > 3)
    {
      cout<<"starting CreateSpline("<<filename<<")"<<endl;
    }

  if(filename.isNull())
    {
      cout<<"Error in "<<__PRETTY_FUNCTION__<<", filename not specified ("<<filename<<")"<<endl;
      exit(2);
    }
  
  if(filename.contains("Rutherford"))
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"filename: "<<filename<<", creating rutherfordian spline"<<endl;
	}

      for(i=0;i<data->NumberOfLevel;i++)
	{
	  TF1* rutherford = new TF1("rutherford","1./TMath::Power(TMath::Sin(x/360.*TMath::Pi()),4.)",0.,180.);
	  //AngularDistribution[i] = new TSpline3("Rutherford",0.,180.,rutherford,1800);
	  AngularDistribution[i] = new TGraph(rutherford);
	  MaxAngularDistribution[i] = AngularDistribution[i]->Eval(ThetaCmMin[i]);
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"created rutherfordian spline (i = "<<i<<"), Maximum = "<<MaxAngularDistribution[i]<<endl;
	    }
	  delete rutherford;
	  CrossSectionIntegral[i] = 1*millibarn;
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"CreateSplines: CrossSectionIntegral["<<i<<"] = "<<CrossSectionIntegral[i]/millibarn<<" mb"<<endl;
	    }
	}
      return;
    }

  if(filename.compareTo("FLAT") == 0)
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"filename: "<<filename<<", creating flat spline"<<endl;
	}

      x = new G4double*[data->NumberOfLevel];
      y = new G4double*[data->NumberOfLevel];
	  
      for(j=0;j<data->NumberOfLevel;j++)
	{
	  MaxAngularDistribution[j] = 1.;
	  x[j] = new G4double[19];
	  y[j] = new G4double[19];

	  for(i=0;i<19;i++)
	    {
	      x[j][i] = i*10.;
	      y[j][i] = 1.;
	    }
	}

      for(j=0;j<data->NumberOfLevel;j++)
	{
	  i=0;
	  while(x[j][i] < ThetaCmMin[j]/degree)
	    {
	      i++;
//	      cout << x[j][i] << y[j][i] << endl;
			if(i > 18)
		{
		  cerr<<"Error in CreateSplines, supposed to use ThetaCmMin = "<<ThetaCmMin[j]/degree<<" degree, but found no x value above this angle"<<endl;
		  exit(2);
		}
	    }
//	  if(data->VerboseLevel > 1)
//	    {
	      cout<<i<<". entry ("<<x[j][i]<<") > ThetaCmMin = "<<ThetaCmMin[j]/degree<<endl;
//	    }

	  CrossSectionIntegral[j] = Integrate(19, x[j],y[j],i)*millibarn;
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"CreateSplines: CrossSectionIntegral["<<j<<"] = "<<CrossSectionIntegral[j]/millibarn<<" mb"<<endl;
	    }
	}

      for(j=0;j<data->NumberOfLevel;j++)
	{
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"creating user defined spline (index = "<<j<<"), Maximum = "<<MaxAngularDistribution[j]<<endl;
	    }
	  
	  //AngularDistribution[j] = new TSpline3(filename.data(),x[j],y[j],entries[j]-skipped_entries[j]);//"e1b1",0,0 => first derivative should be zero at beginning and end of spline
	  AngularDistribution[j] = new TGraph(19,x[j],y[j]);
      
	  if(data->VerboseLevel > 7)
	    {
	      for(i=0;i<19;i++)
		{	
		  //AngularDistribution[j]->GetKnot(i,xtmp,ytmp);
		  AngularDistribution[j]->GetPoint(i,xtmp,ytmp);
		  cout<<i<<": "<<xtmp<<", "<<ytmp<<endl;
		}
	      
	      cout<<"MaxAngularDistribution["<<j<<"] = "<<MaxAngularDistribution[j]<<endl;
	      TFile* f;
	      if(j == 0)
		f = new TFile("tmpfile.root","recreate");
	      else
		f = new TFile("tmpfile.root","update");
	      char tmp[256];
	      sprintf(tmp,"spline%d",j);
	      AngularDistribution[j]->SetName(tmp);
	      AngularDistribution[j]->Write();
	      f->Close();
	    }
	}
      return;
    }

  file.open(filename.data());

  if(!file.is_open())
    {
      cout<<"Sorry, failed to open file "<<filename.data()<<endl;
      exit(1);
    }

  if(file.bad() || file.fail())
    {
      cout<<"Sorry, couldn't open file "<<filename.data()<<endl;
      exit(13);
    }

  if(data->VerboseLevel > 3)
    {
      cout<<"opened file "<<filename<<endl;
    }

  entries = new G4int[data->NumberOfLevel];
  maxentries = 0;
  for(j=0;j<data->NumberOfLevel;j++)
    {
      file>>entries[j];
      if(entries[j] > maxentries)
	maxentries = entries[j];
    }
  file.ignore(1000,'\n');

  if(data->VerboseLevel > 3)
    {
      cout<<"supposed to read ";
      for(j=0;j<data->NumberOfLevel;j++)
	{
	  cout<<entries[j]<<", ";
	}
      cout<<"maxentries = "<<maxentries<<" from "<<filename<<endl;
    }

  x = new G4double*[data->NumberOfLevel];
  y = new G4double*[data->NumberOfLevel];

  G4int* skipped_entries = new G4int[data->NumberOfLevel];

  for(j=0;j<data->NumberOfLevel;j++)
    {
      MaxAngularDistribution[j] = 0;
      x[j] = new G4double[entries[j]+1];//+1 to be able to store additional data point at 180 degree later
      y[j] = new G4double[entries[j]+1];
      skipped_entries[j] = 0;
    }

  for(i=0;i<maxentries;i++)
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"read entry "<<i<<": ";
	}
      for(j=0;j<data->NumberOfLevel;j++)
	{
	  if(i < entries[j])
	    {
	      file>>x[j][i - skipped_entries[j]]>>y[j][i - skipped_entries[j]];
	      if(i - skipped_entries[j] > 0)
		if(x[j][i - skipped_entries[j]] == x[j][i - skipped_entries[j] - 1])
		  {
		    skipped_entries[j]++;
		    if(data->VerboseLevel > 3)
		      {
			cout<<"     skipped      \t";
		      }
		    continue;
		  }
	      y[j][i - skipped_entries[j]] *= sin(x[j][i - skipped_entries[j]]/180.*PI);
	      
	      if(y[j][i - skipped_entries[j]] > MaxAngularDistribution[j])
		MaxAngularDistribution[j] = y[j][i - skipped_entries[j]];
	      if(data->VerboseLevel > 3)
		{
		  cout<<"x = "<<x[j][i - skipped_entries[j]]<<", y = "<<y[j][i - skipped_entries[j]]<<" \t";
		}
	    }
	  else
	    {
	      if(data->VerboseLevel > 3)
		{
		  cout<<"                  \t";
		}
	    }
	}
      cout<<endl;
      file.ignore(1000,'\n');
    }

  file.close();

  if(data->VerboseLevel > 3)
    {
      cout<<"read ";
      for(j=0;j<data->NumberOfLevel;j++)
	{
	  cout<<entries[j]-skipped_entries[j]<<"/"<<skipped_entries[j]<<" ";
	}
      cout<<"entries"<<endl;
    }

  for(j=0;j<data->NumberOfLevel;j++)
    {
      if(x[j][entries[j]-skipped_entries[j]-1] < 180.)
	{
	  x[j][entries[j]-skipped_entries[j]] = 180.;
	  y[j][entries[j]-skipped_entries[j]] = 0.;
	  entries[j]++;
	}

      first=0;
      ifirst=0;
      while(x[j][first] < ThetaCmMin[j]/degree){
	first++;
	if(first > entries[j])
	  {
	    cerr<<"Error in CreateSplines, supposed to use ThetaCmMin = "<<ThetaCmMin[j]/degree<<" degree, but found no x value above this angle"<<endl;
	    exit(2);
	  }
      }
      while(x[j][ifirst] < IntThetaCmMin/degree){
	ifirst++;
	if(ifirst > entries[j])
	  {
	    cerr<<"Error in CreateSplines, supposed to use ThetaCmMin = "<<ThetaCmMin[j]/degree<<" degree, but found no x value above this angle"<<endl;
	    exit(2);
	  }
      }
      
      last = entries[j]-skipped_entries[j]-1;
      ilast = entries[j]-skipped_entries[j]-1;
      while(x[j][last] > ThetaCmMax[j]/degree){
	//cout << "x["<<j<<"]["<<last<<"] " << x[j][last] << " ThetaCmMax["<<j<<"]/degree " << ThetaCmMax[j]/degree << endl;
	last--;
	if(last<1){
	  cerr<<"Error in CreateSplines, supposed to use ThetaCmMax = "<<ThetaCmMax[j]/degree<<" degree, but found no x value below this angle"<<endl;
	  exit(2);
	}
      }
      while(x[j][ilast] > IntThetaCmMax/degree){
	//cout << "x["<<j<<"]["<<ilast<<"]" << x[j][ilast] << " IntThetaCmMax " << IntThetaCmMax/degree << endl;
	ilast--;
	if(ilast<1){
	  cerr<<"Error in CreateSplines, supposed to use ThetaCmMax = "<<ThetaCmMax[j]/degree<<" degree, but found no x value below this angle"<<endl;
	  exit(2);
	}
      }
      if(data->VerboseLevel > -1)
	{
	  cout<<first<<". entry ("<<x[j][first]<<") > ThetaCmMin = "<<ThetaCmMin[j]/degree<<endl;
	  cout<<last<<". entry ("<<x[j][last]<<") < ThetaCmMax = "<<ThetaCmMax[j]/degree<<endl;
	  cout<<ifirst<<". entry ("<<x[j][first]<<") > IntThetaCmMin = "<<IntThetaCmMin/degree<<endl;
	  cout<<ilast<<". entry ("<<x[j][last]<<") < IntThetaCmMax = "<<IntThetaCmMax/degree<<endl;
	}

      //kwimmer
      CrossSectionIntegral[j] = Integrate(last-first, x[j],y[j],first)*millibarn;
      ICrossSectionIntegral[j] += Integrate(ilast-ifirst, x[j],y[j],ifirst)*millibarn;
      if(data->VerboseLevel > -3)
	{
	  cout<<"CreateSplines: CrossSectionIntegral["<<j<<"] = "<<CrossSectionIntegral[j]/millibarn<<" mb"<<endl;
	  cout << "ICrossSectionIntegral["<<j<<"] = "<<ICrossSectionIntegral[j]/millibarn<<" mb"<<endl;
	}
    }

  for(j=0;j<data->NumberOfLevel;j++)
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"creating user defined spline (index = "<<j<<"), Maximum = "<<MaxAngularDistribution[j]<<endl;
	}

      //AngularDistribution[j] = new TSpline3(filename.data(),x[j],y[j],entries[j]-skipped_entries[j]);//"e1b1",0,0 => first derivative should be zero at beginning and end of spline
      //DownScaling 
      //cout << "j " << j << " DownScale " << data->DownScale << endl;
      if(j==0&&data->DownScale>0&&data->Coulex){
	//cout << "down scaling with " << data->DownScale << endl;
	for(i=0; i<entries[0]-skipped_entries[0];i++){
	  // cout << y[0][i] << " -> " ;
	  y[0][i]/=(double)data->DownScale;
	  //cout << y[0][i] << endl;
	}
      }
      AngularDistribution[j] = new TGraph(entries[j]-skipped_entries[j],x[j],y[j]);
      
      if(data->VerboseLevel > 3)
	{
	  cout << "-----------------------------------------------------" << endl;
	  for(i=0;i<maxentries;i++)
	    {	
	      if(i < entries[j] - skipped_entries[j])
		{
		  //AngularDistribution[j]->GetKnot(i,xtmp,ytmp);
		  AngularDistribution[j]->GetPoint(i,xtmp,ytmp);
		  cout<<i<<": "<<xtmp<<", "<<ytmp<<endl;
		}
	    }

	  cout<<"MaxAngularDistribution["<<j<<"] = "<<MaxAngularDistribution[j]<<endl;
	  /*
	  TFile* f;
	  if(j == 0)
	    f = new TFile(stmp.data(),"recreate");
	  else
	    f = new TFile(stmp.data(),"update");
	  char tmp[256];
	  sprintf(tmp,"spline%d",j);
	  AngularDistribution[j]->SetName(tmp);
	  AngularDistribution[j]->Write();
	  f->Close();
	  cout << "-----------------------------------------------------" << endl;
	  */
	}
    }

  delete[] entries;
  delete[] skipped_entries;

  //possible without destroying TGraph ???
  for(j=0;j<data->NumberOfLevel;j++)
    {
      delete[] x[j];
      delete[] y[j];
    }
  delete[] x;
  delete[] y;

  return;
}

void PrimaryGeneratorAction::CreateOpticalPotential(G4String filename)
{
  G4int i, j;
  G4int entries = 0;
  ifstream file;
  G4double* x;
  G4double** y;

  G4String stmp = data->SimulationDirectory+"/OpticalPotential.root";

  if(data->VerboseLevel > 3)
    {
      cout<<"starting CreateOpticalPotential("<<filename<<")"<<endl;
    }

  if(filename.isNull())
    {
      cout<<"Error in "<<__PRETTY_FUNCTION__<<", filename not specified ("<<filename<<")"<<endl;
      exit(2);
    }

  if(filename.contains("Rutherford"))
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"filename: "<<filename<<", creating rutherfordian spline"<<endl;
	}

      char tmp[256];

      for(i=0;i<data->TargetMaterial->NumberOfElements();i++)
	{
	  sprintf(tmp,"1./TMath::Power(TMath::Sin(x/360.*TMath::Pi()),4.)*%g/%g",RutherfordFactor[i]/(pow(MeV,2)*millibarn), pow(E_cm[i]/MeV,2));
	  TF1* rutherford = new TF1("rutherford",tmp,0.,180.);
	  //OpticalPotential[i] = new TSpline3("Rutherford",0.,180.,rutherford,1800);
	  OpticalPotential[i] = new TGraph(rutherford);
	  MaxOpticalPotential[i] = OpticalPotential[i]->Eval(ThetaCmMin[i+data->NumberOfLevel]);
	  //CrossSectionIntegral[i+data->NumberOfLevel] = RutherfordFactor[i]/(pow(MeV,2)*millibarn)/pow(E_cm[i]/MeV,2)*2./pow(tan(ThetaCmMin[i+data->NumberOfLevel]/360.*PI),2);//*millibarn;
	  CrossSectionIntegral[i+data->NumberOfLevel] = 1*millibarn;
	  if(data->VerboseLevel > 1)
	    {
	      cout<<"created rutherfordian spline (i = "<<i<<"), Maximum = "<<MaxOpticalPotential[i]<<endl;
	      cout<<"CreateOpticalPotential: CrossSectionIntegral["<<i+data->NumberOfLevel<<"] = "<<CrossSectionIntegral[i+data->NumberOfLevel]/millibarn<<" mb = "<<RutherfordFactor[i]/(pow(MeV,2)*millibarn)<<"*"<<1/pow(E_cm[i]/MeV,2)<<"*"<<2./pow(tan(ThetaCmMin[i+data->NumberOfLevel]/360.*PI),2)<<endl;
	    }
	  delete rutherford;
	}
      return;
    }

  file.open(filename.data());

  if(!file.is_open())
    {
      cout<<"Sorry, failed to open file "<<filename.data()<<endl;
      exit(1);
    }

  if(file.bad() || file.fail())
    {
      cout<<"Sorry, couldn't open file "<<filename.data()<<endl;
      exit(13);
    }

  if(data->VerboseLevel > 3)
    {
      cout<<"opened file "<<filename<<endl;
    }

  file.ignore(1000,'\n');
  file.ignore(1000,'\n');
  file>>data->OpticalFactor>>data->OpticalExponent;
  file.ignore(1000,'\n');

  if(data->VerboseLevel > 3)
    {
      cout<<filename<<": Factor = "<<data->OpticalFactor<<", Exponent = "<<data->OpticalExponent<<endl;
    }

  file>>entries;
  file.ignore(1000,'\n');

  if(data->VerboseLevel > 3){
    cout<<"supposed to read "<<entries<<" entries from "<<filename<<endl;
    cout<<"number of columns "<<data->TargetMaterial->NumberOfElements()<<endl;
  }
  x = new G4double[entries+1];
  y = new G4double*[data->TargetMaterial->NumberOfElements()];
  for(j=0;j<data->TargetMaterial->NumberOfElements();j++){
    y[j] = new G4double[entries+1];
    MaxOpticalPotential[j] = 0;
  }
  for(i=0;i<entries;i++){
    file>>x[i];
    for(j=0;j<data->TargetMaterial->NumberOfElements();j++){
      file>>y[j][i];
      y[j][i] *= sin(x[i]/180.*PI); // commented out 15.9	  
      if(y[j][i] > MaxOpticalPotential[j])
	MaxOpticalPotential[j] = y[j][i];
    }
    if(data->VerboseLevel > 3){
      cout<<"read entry "<<i<<"/"<<entries<<": x = "<<x[i];
      for(j=0;j<data->TargetMaterial->NumberOfElements();j++)
	cout<<", y["<<j<<"] = "<<y[j][i];
      cout<<endl;
    }
    file.ignore(1000,'\n');
  }
  
  file.close();
  
  if(x[entries-1] < 180.){
    x[entries] = 180.;
    for(j=0;j<data->TargetMaterial->NumberOfElements();j++)
      y[j][entries] = 0.;
    entries++;
  }
  
  for(j=0;j<data->TargetMaterial->NumberOfElements();j++){
    i=0;
    while(x[i] < ThetaCmMin[j+data->NumberOfLevel]/degree){
      i++;
      if(i > entries){
	cerr<<"Error in CreateOpticalPotential, supposed to use ThetaCmMin = "<<ThetaCmMin[j+data->NumberOfLevel]/degree<<" degree, but found no x value above this angle"<<endl;
	exit(2);
      }
    }
    if(data->VerboseLevel > 1){
      cout<<i<<". entry ("<<x[i]<<") > ThetaCmMin = "<<ThetaCmMin[j+data->NumberOfLevel]/degree<<endl;
    }
    
    CrossSectionIntegral[j+data->NumberOfLevel] = Integrate(entries, x, y[j],i)*millibarn;
    
    if(data->VerboseLevel > 3){
      cout<<"CreateOpticalPotential: CrossSectionIntegral["<<j+data->NumberOfLevel<<"] = "<<CrossSectionIntegral[j+data->NumberOfLevel]/millibarn<<" mb"<<endl;
      cout<<"creating user defined spline (index = "<<j<<"), Maximum = "<<MaxOpticalPotential[j]<<endl;
    }
    
    //OpticalPotential[j] = new TSpline3(filename.data(),x,y[j],entries);
    OpticalPotential[j] = new TGraph(entries,x,y[j]);
    /* doesn't work
       if(data->VerboseLevel > 3){
       cout<<"MaxOpticalPotential["<<j<<"] = "<<MaxOpticalPotential[j]<<endl;
       TFile* f;
       if(j == 0)
       f = new TFile(stmp.data(),"recreate");
       else
       f = new TFile(stmp.data(),"update");
       //f.cd();
       OpticalPotential[j]->Write(Form("OpticalPotential%d",j),TObject::kOverwrite);
       f->Close();
       }
    */
  }
  

  //possible without destroying TGraph ???
  //delete[] x;
  //delete[] y;

  return ;
}

G4double PrimaryGeneratorAction::Integrate(G4int entries, G4double* x, G4double* y, G4int first)
{
  if(data->VerboseLevel > 1)
    {
      cout<<"starting integration of "<<entries<<" entries, starting with "<<first<<". entry"<<endl;
    }

  G4double Integral = 0;

//  for(int i=first+1;i<entries+first;i++)
  for(int i=first+1;i<entries;i++)
    {
      if(x[i] < 0. || 180. < x[i])
	{
	  cout<<"x-input "<<x[i]<<" is out of allowed range (0-180 deg.)"<<endl;
	  return -1.;
	}
      //Integral += (y[i]*sin(x[i]/180.*PI)+y[i-1]*sin(x[i-1]/180.*PI))/2.*(x[i]-x[i-1]);//original vinzenz
      Integral += (y[i]+y[i-1])/2.*(x[i]-x[i-1])/180.*PI; //modified 25.03.10


      //cout <<"x["<<i<<"] = " << x[i] << " Integral["<<i<<"]: " << Integral <<" = "<< endl;
      // cout << "(y[i]*sin(x[i]/180.*PI)\t +y[i-1]*sin(x[i-1]/180.*PI))\t /2.*(x[i]-x[i-1]) = " << endl;
      //cout << y[i]<<"*"<<sin(x[i]/180.*PI)<<"\t +" << y[i-1]<<"*"<<sin(x[i-1]/180.*PI) << "\t " << 0.5*(x[i]-x[i-1]) << endl;
      //cout << "x[i] " << x[i] << "x[i-1] " << x[i-1] << endl;
    }

  if(data->VerboseLevel > 1)
    {
      cout<<"integration of "<<entries<<" entries from "<<x[first]<<" to "<<x[first+entries-1]<<" yields "<<2.*PI*Integral<<endl;
    }

  return 2.*PI*Integral;
}

void PrimaryGeneratorAction::BeamParticles()
{
  if(!data->Efficiency && !data->AlphaSource && !data->QuadAlphaSource && !data->BetaEnergy)
    {
      data->BeamParticles = data->EventNumber/ScatteringProbability[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()-1];
      if(data->Coulex)
	data->BeamParticles = data->EventNumber/ScatteringProbability[data->NumberOfLevel];
      //if(data->VerboseLevel > 1)
      //{
      if(!data->Coulex)
	cout<<"BeamParticles = EventNumber/ScatteringProbability = "<<data->EventNumber<<"/"<<ScatteringProbability[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()-1]<<" = "<<data->BeamParticles<<endl;
      else
	cout<<"BeamParticles = EventNumber/ScatteringProbability = "<<data->EventNumber<<"/"<<ScatteringProbability[data->NumberOfLevel]<<" = "<<data->BeamParticles<<endl;
	
      //}
    }

  return;
}

void PrimaryGeneratorAction::SetThetaLab(double ThetaLab[])
{
  if(data->ExperimentType == 0)
    {
      ThetaLab[0] = atan(data->settings->SSDActInnerRad/data->settings->SSDZPosition)*radian;
      ThetaLab[1] = atan(data->settings->SSDActOuterRad/data->settings->SSDZPosition)*radian;
      ThetaLab[2] = ThetaLab[0];
      ThetaLab[3] = ThetaLab[1];
      if(data->VerboseLevel > 3)
	{
	  cout<<"calculated ThetaLab from atan of "<<(data->settings->SSDActInnerRad/data->settings->SSDZPosition)<<" and "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))<<" ("<<data->settings->ForwardBarrelZ<<" - "<<data->settings->BarrelActLength/2.<<")"<<endl;
	}
    }
  else if(data->ExperimentType == 4)
    {
      ThetaLab[0] = atan(data->settings->SSDActInnerRad/data->settings->SSDZPosition)*radian;
      ThetaLab[1] = acos((data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.)/sqrt(pow(data->settings->BarrelDistance,2)+pow(data->settings->BarrelActLength/2.,2)+pow(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.,2))); //including edges
      //atan(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))*radian;
      ThetaLab[2] = PI - ThetaLab[1];
      ThetaLab[3] = PI - ThetaLab[0];
		if(data->Coulex){
			ThetaLab[2] = ThetaLab[0];
			ThetaLab[3] = ThetaLab[1];
      }
      if(data->VerboseLevel > 3)
	{
	  cout<<"calculated ThetaLab from atan of "<<(data->settings->SSDActInnerRad/data->settings->SSDZPosition)<<" and "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))<<" ("<<data->settings->ForwardBarrelZ<<" - "<<data->settings->BarrelActLength/2.<<")"<<endl;
	}
    }
  else if(data->ExperimentType == 1)
    {
      ThetaLab[0] = PI - atan(data->settings->SSDActOuterRad/data->settings->SSDZPosition)*radian;
      ThetaLab[1] = PI - atan(data->settings->SSDActInnerRad/data->settings->SSDZPosition)*radian;
      ThetaLab[2] = ThetaLab[0];
      ThetaLab[3] = ThetaLab[1];
      if(data->VerboseLevel > 3)
	{
	  cout<<"calculated ThetaLab from atan of "<<(data->settings->SSDActOuterRad/data->settings->SSDZPosition)<<" and "<<(data->settings->SSDActInnerRad/data->settings->SSDZPosition)<<endl;
	}
    }
  else if(data->ExperimentType == 2)
    {
      ThetaLab[0] = atan(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ+data->settings->BarrelActLength/2.))*radian;
      ThetaLab[1] = atan(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))*radian;
      ThetaLab[2] = PI - ThetaLab[1];
      ThetaLab[3] = PI - ThetaLab[0];
      if(data->VerboseLevel > 3)
	{
	  cout<<"calculated ThetaLab from atan of "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ+data->settings->BarrelActLength/2.))<<" and "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))<<" ("<<data->settings->ForwardBarrelZ<<" +- "<<data->settings->BarrelActLength/2.<<")"<<endl;
	}
    }
  else if(data->ExperimentType == 3)
    {
      ThetaLab[0] = atan(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ+data->settings->BarrelActLength/2.))*radian;
      ThetaLab[1] = atan(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))*radian;
      ThetaLab[2] = PI - ThetaLab[1];
      ThetaLab[3] = PI - atan(data->settings->SSDActInnerRad/data->settings->SSDZPosition)*radian;
      if(data->VerboseLevel > 3)
	{
	  cout<<"calculated ThetaLab from atan of "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ+data->settings->BarrelActLength/2.))<<" and "<<(data->settings->BarrelDistance/(data->settings->ForwardBarrelZ-data->settings->BarrelActLength/2.))<<" and "<<(data->settings->SSDActInnerRad/data->settings->SSDZPosition)<<" ("<<data->settings->ForwardBarrelZ<<" +- "<<data->settings->BarrelActLength/2.<<")"<<endl;
	}
    }
  else
    {
      cerr<<"Experiment type "<<data->ExperimentType<<" is not implemented yet, exiting"<<endl;
      exit(3);
    }

  return;
}

void PrimaryGeneratorAction::MinimalThetaCm(double ThetaLab[], double energy[])
{
  int i;
  double tmp[4]={0.,0.,0.,0.};
  if(data->Coulex)
    ThetaCmMin = new double[data->NumberOfLevel];
  else
    ThetaCmMin = new double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];

  //for the different levels populated in the transfer reaction
  if(!data->Coulex){
    for(i=0;i<data->NumberOfLevel;i++){
      if(data->ThetaCmMin > 0.*degree)
	ThetaCmMin[i] = 0.*degree;//these are transfer reactions
      else{
	if(data->Projectile->Mass() > data->Target->Mass()){
	  ThetaCmMin[i] = kinematics->CMAngle(ThetaLab[3]/degree,data->Projectile,data->Target,data->Recoil,data->Ejectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;
	}
	else{
	  ThetaCmMin[i] = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Recoil,data->Ejectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;
	}
	if(data->VerboseLevel > 3){
	  cout<<"Original: Transfer reaction of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->Target->A()<<data->Target->Name()<<", E_ex = "<<data->LevelEnergy[i]/MeV<<" MeV: theta_cm_min = "<<ThetaCmMin[i]/degree<<" degree (ThetaLab = "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	if(ThetaCmMin[i] > 1.*degree)
	  ThetaCmMin[i] -= 1*degree;
	else
	  ThetaCmMin[i] = 0.*degree;
      }
    }
    //for the elastic scattering
    for(i=data->NumberOfLevel;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++){
      if(data->ThetaCmMin > 0.*degree)
	ThetaCmMin[i] = data->ThetaCmMin;
      else{
	ThetaCmMin[i] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->Projectile,2,energy,0.,0)*degree;
	tmp[0]           = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->Projectile,2,energy,0.,1)*degree;
	
	if(tmp[0] < ThetaCmMin[i] && tmp[0] >= 0)//tmp < 0 means there is no solution, i.e. the particle never hits the detector
	  ThetaCmMin[i] = tmp[0];
	
	if(ThetaCmMin[i] < 0){//means there is no solution, i.e. the particle never hits the detector
	  
	  if(tmp[0] >= 0)
	    ThetaCmMin[i] = tmp[0];
	  else
	    cerr<<"Error, both solutions are not valid ("<<ThetaCmMin[i]/degree<<", "<<tmp[0]/degree<<")"<<endl;
	}
	
	if(data->VerboseLevel > 3){
	  cout<<"Original: Elastic Scattering of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->A()<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->Name()<<": theta_cm_min = "<<ThetaCmMin[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	
	if(ThetaCmMin[i] > 1.*degree)
	  ThetaCmMin[i] -= 1*degree;
	else
	  ThetaCmMin[i] = 0.*degree;
      }
    }
  } 
  else{//coulex
    for(i=0;i<data->NumberOfLevel;i++){
      if(data->ThetaCmMin > 0.*degree){
	ThetaCmMin[i] = data->ThetaCmMin*degree;
	IntThetaCmMin = data->ThetaCmMin*degree;
      }
      else{
	//cout << "-------------------------------- " << endl;
	//cout << "start thmin " << endl;
	tmp[0] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,0)*degree;
	//cout << ThetaLab[0]/degree << " and code 0 " << tmp[0]/degree << endl; 
	tmp[1] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,1)*degree;
	//cout << ThetaLab[0]/degree << " and code 1 " << tmp[1]/degree << endl; 
	tmp[2] = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,0)*degree;
	//cout << ThetaLab[1]/degree << " and code 0 " << tmp[2]/degree << endl; 
	tmp[3] = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,1)*degree; //level energy always 0 for coulex
	//cout << ThetaLab[1]/degree << " and code 1 " << tmp[3]/degree << endl; 
	for(int k=0;k<4;k++){
	  if(tmp[k]<0)
	    tmp[k] = 1000;
	}
	ThetaCmMin[i] = min(tmp[0],tmp[1]);
	ThetaCmMin[i] = min(ThetaCmMin[i],tmp[2]);
	ThetaCmMin[i] = min(ThetaCmMin[i],tmp[3]);
	IntThetaCmMin = min(tmp[1],tmp[3]); //recoil angles!
	cout << " __________________________" << endl;
	cout << " angles " << tmp[0]/degree << ", " << tmp[1]/degree << ", " << tmp[2]/degree << ", " << tmp[3]/degree << "-> min " << ThetaCmMin[i]/degree<<endl;
	/*
	ThetaCmMin[i] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,0)*degree;

	tmp           = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;
	cout << " __________________________" << endl;
	cout << " 0 0 " << kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,0) << endl;
	cout << " 0 1 " << kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,1) << endl;
	cout << " 1 0 " << kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,0) << endl;
	cout << " 1 1 " << kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,1) << endl;
	cout << " __________________________" << endl;
	if(tmp < ThetaCmMin[i] && tmp >= 0)//tmp < 0 means there is no solution, i.e. the particle never hits the detector
	  ThetaCmMin[i] = tmp;
	
	if(ThetaCmMin[i] < 0){//means there is no solution, i.e. the particle never hits the detector
	  
	  if(tmp >= 0)
	    ThetaCmMin[i] = tmp;
	  else
	    cerr<<"Error, both solutions are not valid ("<<ThetaCmMin[i]/degree<<", "<<tmp/degree<<")"<<endl;
	}
	*/
	if(data->VerboseLevel > 3){
	  cout << "output " << endl;
	  cout<<"Original: Coulex of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->Target->A()<<data->Target->Name()<<"eexc " <<data->LevelEnergy[i]/MeV<<": theta_cm_min = "<<ThetaCmMin[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	
	if(ThetaCmMin[i] > 1.*degree)
	  ThetaCmMin[i] -= 1*degree;
	else
	  ThetaCmMin[i] = 0.*degree;
      }
    }
  }
  return;
}
void PrimaryGeneratorAction::MaximalThetaCm(double ThetaLab[], double energy[])
{
  int i;
  double tmp[4]={0.,0.,0.,0.};
  if(data->Coulex)
    ThetaCmMax = new double[data->NumberOfLevel];
  else
    ThetaCmMax = new double[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()];

  //for the different levels populated in the transfer reaction
  if(!data->Coulex){
    for(i=0;i<data->NumberOfLevel;i++){
      if(data->ThetaCmMax < 180.*degree)
	ThetaCmMax[i] = 180.*degree;//these are transfer reactions
      else{
	if(data->Projectile->Mass() > data->Target->Mass()){
	  ThetaCmMax[i] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Recoil,data->Ejectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;
	}
	else{
	  ThetaCmMax[i] = kinematics->CMAngle(ThetaLab[2]/degree,data->Projectile,data->Target,data->Recoil,data->Ejectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;
	}
	if(data->VerboseLevel > 3){
	  cout<<"Original: Transfer reaction of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->Target->A()<<data->Target->Name()<<", E_ex = "<<data->LevelEnergy[i]/MeV<<" MeV: theta_cm_max = "<<ThetaCmMax[i]/degree<<" degree (ThetaLab = "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	if(ThetaCmMax[i] > 1.*degree)
	  ThetaCmMax[i] += 1*degree;
	else
	  ThetaCmMax[i] = 0.*degree;
      }
    }
    //for the elastic scattering
    for(i=data->NumberOfLevel;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++){
      if(data->ThetaCmMax < 180.*degree)
	ThetaCmMax[i] = data->ThetaCmMax;
      else{
	ThetaCmMax[i] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->Projectile,2,energy,0.,1)*degree;
	tmp[0]           = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->TargetMaterial->GetElement(i-data->NumberOfLevel),data->Projectile,2,energy,0.,0)*degree;
	
	if(tmp[0] > ThetaCmMax[i] && tmp[0] >= 0)//tmp < 0 means there is no solution, i.e. the particle never hits the detector
	  ThetaCmMax[i] = tmp[0];
	
	if(ThetaCmMax[i] < 0){//means there is no solution, i.e. the particle never hits the detector
	  
	  if(tmp[0] >= 0)
	    ThetaCmMax[i] = tmp[0];
	  else
	    cerr<<"Error, both solutions are not valid ("<<ThetaCmMax[i]/degree<<", "<<tmp[0]/degree<<")"<<endl;
	}
	
	if(data->VerboseLevel > 3){
	  cout<<"Original: Elastic Scattering of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->A()<<data->TargetMaterial->GetElement(i-data->NumberOfLevel)->Name()<<": theta_cm_max = "<<ThetaCmMax[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	
	if(ThetaCmMax[i] > 1.*degree)
	  ThetaCmMax[i] += 1*degree;
	else
	  ThetaCmMax[i] = 0.*degree;
      }
    }
  } 
  else{//coulex
    for(i=0;i<data->NumberOfLevel;i++){
      if(data->ThetaCmMax < 180){
	cout << data->ThetaCmMax*degree << endl;
	ThetaCmMax[i] = data->ThetaCmMax*degree;
	IntThetaCmMax = data->ThetaCmMax*degree;
      }
      else{
	//cout << "-------------------------------- " << endl;
	//cout << "start thmax " << endl;
	
	tmp[0] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,0)*degree;
	//cout << ThetaLab[0]/degree << " and code 0 " << tmp[0]/degree << endl; 
	tmp[1] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,1)*degree;
	//cout << ThetaLab[0]/degree << " and code 1 " << tmp[1]/degree << endl; 
	tmp[2] = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,0)*degree;
	//cout << ThetaLab[1]/degree << " and code 0 " << tmp[2]/degree << endl; 
	tmp[3] = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,0.0/MeV,1)*degree; // level energy always 0 for coulex
	//cout << ThetaLab[1]/degree << " and code 1 " << tmp[3]/degree << endl; 
	for(int k=0;k<4;k++){
	  if(tmp[k]<0)
	    tmp[k] = -1000;
	}
	ThetaCmMax[i] = max(tmp[0],tmp[1]);
	ThetaCmMax[i] = max(ThetaCmMax[i],tmp[2]);
	ThetaCmMax[i] = max(ThetaCmMax[i],tmp[3]);
	IntThetaCmMax = max(tmp[1],tmp[3]);
	cout << " __________________________" << endl;
	cout << " angles " << tmp[0]/degree << ", " << tmp[1]/degree << ", " << tmp[2]/degree << ", " << tmp[3]/degree << "-> max " << ThetaCmMax[i]/degree<<endl;

	/*
	ThetaCmMax[i] = kinematics->CMAngle(ThetaLab[0]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,1)*degree;

	tmp           = kinematics->CMAngle(ThetaLab[1]/degree,data->Projectile,data->Target,data->Target,data->Projectile,2,energy,data->LevelEnergy[i]/MeV,0)*degree;
	
	//cout << "ThetaLab[0] = " << ThetaLab[0]/degree << " in cm " << ThetaCmMax[i]/degree << endl;
	//cout << "ThetaLab[1] = " << ThetaLab[1]/degree << " in cm " << tmp/degree << endl;

	if(tmp > ThetaCmMax[i] && tmp >= 0)//tmp < 0 means there is no solution, i.e. the particle never hits the detector
	  ThetaCmMax[i] = tmp;
	
	if(ThetaCmMax[i] < 0){//means there is no solution, i.e. the particle never hits the detector
	  
	  if(tmp >= 0)
	    ThetaCmMax[i] = tmp;
	  else{
	    cerr<<"Error, both solutions are not valid ("<<ThetaCmMax[i]/degree<<", "<<tmp/degree<<")"<<endl;
	    ThetaCmMax[i] = 180.*degree;
	  }
	}
	*/
	if(data->VerboseLevel > 3){
	  cout << "output " << endl;
	  cout<<"Original: Coulex of "<<data->Projectile->A()<<data->Projectile->Name()<<" on "<<data->Target->A()<<data->Target->Name()<<"eexc " <<data->LevelEnergy[i]/MeV<<": theta_cm_max = "<<ThetaCmMax[i]/degree<<" degree (ThetaLab = "<<ThetaLab[0]/degree<<" or "<<ThetaLab[1]/degree<<" or "<<ThetaLab[3]/degree<<" degree)"<<endl;
	}
	
	if(ThetaCmMax[i] > 1.*degree)
	  ThetaCmMax[i] += 1*degree;
	else
	  ThetaCmMax[i] = 0.*degree;
      }
    }
  }
  return;
}
int PrimaryGeneratorAction::ReactionType()
{
  G4double tmp;
  int reaction;
  if(data->IsotropeParticleDistribution)
    {
      reaction = 0;
      data->Transfer = 1;
    }

  else if(data->TransferCrossSection >= 0)
    {
      //decide whether to simulate transfer reaction or elastic scattering
      tmp = RandFlat::shoot(0.,1.);
		if(tmp < data->TransferCrossSection) //transfer reaction
//		if (1)
		{
	  data->Transfer = 1;
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"shooting from 0 to "<<ScatteringProbability[data->NumberOfLevel-1]<<endl;
	    }
	  tmp = RandFlat::shoot(0.,ScatteringProbability[data->NumberOfLevel-1]);
	      
	  for(reaction=0;tmp > ScatteringProbability[reaction];reaction++)
	    {
	      if(reaction >= data->NumberOfLevel)
		{
		  cerr<<"Error in PrimaryGeneratorAction: couldn't find reaction channel for tmp = "<<tmp<<endl;
		  exit(1);
		}
	      //nothing else, after this loop reaction is the reaction channel to be used (from 0 to NumberOfLevel-1 or NumberOfLevel+NumberOfElements-1)
	    }
	}
      else//elastic scattering => reaction based on cross sections for elastic scattering on the different elements in the target
	{
	  data->Transfer = 0;
	  if(data->TargetMaterial->NumberOfElements() == 1){ //easy, only one element in target => chose element 0
	    reaction = data->NumberOfLevel;
	  }
	  else
	    {
	      if(data->noCarbon ==0){
		tmp = RandFlat::shoot(0.,1.);
		for(reaction=data->NumberOfLevel;reaction<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();reaction++)
		  {
		    if(tmp < (part_sum[reaction-data->NumberOfLevel]/sum))
		      break;
		  }
	      }
	      else 
		reaction = data->NumberOfLevel +1;
	    }
	}
    }//if(data->TransferCrossSection >= 0)
  else
    {
      if(data->VerboseLevel > 3)
	{
	  cout<<"shooting from 0 to "<<ScatteringProbability[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()-1]<<endl;
	}
      tmp = RandFlat::shoot(0.,ScatteringProbability[data->NumberOfLevel+data->TargetMaterial->NumberOfElements()-1]);

      for(reaction=0;tmp > ScatteringProbability[reaction];reaction++)
	{
	  if(reaction >= data->NumberOfLevel+data->TargetMaterial->NumberOfElements())
	    {
	      cerr<<"Error in PrimaryGeneratoraction: couldn't find reaction channel for tmp = "<<tmp<<endl;
	      exit(1);
	    }
	  //nothing else, after this loop reaction is the reaction channel to be used (from 0 to NumberOfLevel-1 or NumberOfLevel+NumberOfElements-1)
	}

      //set transfer flag
      if(reaction < data->NumberOfLevel)
	{
	  data->Transfer = 1;
	}
      else
	{
	  data->Transfer = 0;
	}
      if(data->VerboseLevel > 3)
	{
	  cout<<"chose reaction channel "<<reaction<<", tmp = "<<tmp<<", ";
	  for(int i=0;i<data->NumberOfLevel+data->TargetMaterial->NumberOfElements();i++)
	    {
	      cout<<ScatteringProbability[i]<<", ";
	    }
	  cout<<"TransferCrossSection = "<<data->TransferCrossSection<<", Transfer = "<<data->Transfer<<", NumberOfLevel = "<<data->NumberOfLevel<<endl;
	}
    }
  //cout << "data->Transfer " << data->Transfer << endl;
  //cout << reaction << endl;
  return reaction;
}

double PrimaryGeneratorAction::ExcitationEnergy(int reaction)
{
  double eex;


  if(data->Transfer == 1)
    {
//     	cout << "H: " << data->LevelWidth[reaction] << ",	" << data->LevelEnergy[reaction] << endl;

		if (data->LevelWidth[reaction]>0){
			eex = gRandom->BreitWigner(data->LevelEnergy[reaction],data->LevelWidth[reaction]);
			while ((20.*data->LevelWidth[reaction]<abs(eex-data->LevelEnergy[reaction])) || (eex<0)){
				eex = gRandom->BreitWigner(data->LevelEnergy[reaction],data->LevelWidth[reaction]);
			}
			data->EexEnergy = eex;
		}
      else{
			eex = data->LevelEnergy[reaction];
      	data->EexEnergy = eex;
//			cout << eex << endl;
    	}
	 }
  else
    {
      eex = 0.;
      
      data->EexEnergy = -1;
    }
	

  return eex;
}

void PrimaryGeneratorAction::ThetaCm(int reaction, double eex)
{
  double tmp;

  if(data->IsotropeParticleDistribution)
    {
      if(data->VerboseLevel > 4)
	{
	  cout<<"shooting theta_cm 'flat' from 0 degree to 180 degree (eex = "<<eex<<")"<<endl;
	}
      data->Theta_cm = RandFlat::shoot(-1., 1.);
      data->Theta_cm = acos(data->Theta_cm)*radian;
    }
  //use rutherford angular distribution if no ang. distr. file provided (Factor = 0) or if elastic scattering and no optical potential scattering is provided
  else if(data->Factor[reaction] == 0 || (data->Transfer == 0 && !data->OpticalPotential)) 
    {
      if(data->VerboseLevel > 4)
	{
	  cout<<"shooting theta_cm from "<<ThetaCmMin[reaction]/radian<<" radian to "<<M_PI<<" radian (eex = "<<eex<<", reaction = "<<reaction<<", ruth.)"<<endl;
	}
      data->Theta_cm = RandFlat::shoot(-2./3.*(cos(ThetaCmMin[reaction]/radian/2.)/pow(sin(ThetaCmMin[reaction]/radian/2.),3.)+2.*cos(ThetaCmMin[reaction]/radian/2.)/sin(ThetaCmMin[reaction]/radian/2.)),0.);
      tmp = pow(sqrt(1.+pow(0.75*data->Theta_cm,2.))-0.75*data->Theta_cm,1./3.);
      if(data->VerboseLevel > 4)
	{
	  cout<<"theta_cm = RandFlat::shoot("<<-2./3.*(cos(ThetaCmMin[reaction]/radian/2.)/pow(sin(ThetaCmMin[reaction]/radian/2.),3.)+2.*cos(ThetaCmMin[reaction]/radian/2.)/sin(ThetaCmMin[reaction]/radian/2.))<<",0.) = "<<data->Theta_cm<<", tmp = pow(sqrt(1.+"<<pow(0.75*data->Theta_cm,2.)<<")-"<<0.75*data->Theta_cm<<",1./3.) = "<<tmp<<endl;
	}
      data->Theta_cm = -2.*atan((tmp/(1.-pow(tmp,2.)))/radian);
      if(data->Theta_cm < 0.*degree || 180.*degree < data->Theta_cm)
	{
	  cout<<"Error in GeneratePrimaryAction, reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" degree not possible (theta_cm = -2.*atan("<<tmp<<"/"<<1.-pow(tmp,2.)<<"))"<<endl;
	  exit(15);
	}
      if(data->VerboseLevel > 4)
	{
	  file<<data->Theta_cm/degree<<" "<<tmp<<" "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<endl;
	  cout<<"theta_cm = "<<data->Theta_cm/degree<<" degree, tmp = "<<tmp<<", AngDist. = "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<endl;
	}
    }
  //used provided angular distributions (transfer or elastic scattering)
  else
    {
      //transfer reaction with user provided angular distribution
      if(data->Transfer == 1)
	{
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"shooting theta_cm from "<<ThetaCmMin[reaction]/degree<<" radian to "<<ThetaCmMax[reaction]/degree<<" radian (eex = "<<eex<<", reaction = "<<reaction<<")"<<endl;
	    }
	  if(data->VerboseLevel > 3)
	    {
	      cout<<"data->Transfer = "<<data->Transfer<<", data->Factor["<<reaction<<"] = "<<data->Factor[reaction]<<", data->OpticalPotential = "<<data->OpticalPotential<<": starting loop do{}while(tmp > AngularDistribution["<<reaction<<"]->Eval(theta_cm/degree))"<<endl;
	    }
	  do {
	    //new version from theta_cm_min to 180 degree, according to NR
	    //check whether exponential approximation is used
	    //kwimmer 12.3.10 from thcmmin to thcmmax
	    if(data->Exponent[reaction] > 0.)
	      {
		data->Theta_cm = (-log(RandFlat::shoot(exp(-data->Exponent[reaction]*ThetaCmMin[reaction]/degree),exp(-data->Exponent[reaction]*ThetaCmMax[reaction]/degree)))/data->Exponent[reaction])*degree;
		if(data->Theta_cm < ThetaCmMin[reaction] || ThetaCmMax[reaction] < data->Theta_cm)
		  {
		    cout<<"Error in GeneratePrimaryAction, reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" not possible (theta_cm = ((log("<<data->Factor[reaction]<<"/"<<data->Exponent[reaction]<<") - log(RandFlat::shoot("<<data->Factor[reaction]/data->Exponent[reaction]<<","<<data->Factor[reaction]/data->Exponent[reaction]*exp(-180.*data->Exponent[reaction])<<")))/"<<data->Exponent[reaction]<<")*degree"<<endl;
		    exit(15);
		  }
		tmp = RandFlat::shoot(0.,data->Factor[reaction]*exp(-data->Exponent[reaction]*data->Theta_cm/degree));
		if(data->VerboseLevel > 4)
		  {
		    cout<<"GeneratePrimaryAction: reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" degree (RandFlat::shoot("<<exp(-data->Exponent[reaction]*ThetaCmMin[reaction]/degree)<<","<<exp(-180.*data->Exponent[reaction])<<")), tmp = "<<tmp<<" (0-"<<data->Factor[reaction]*exp(-data->Exponent[reaction]*data->Theta_cm/degree)<<"), f(theta_cm) = "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<endl;
		    cout<<"theta_cm = "<<data->Theta_cm/degree<<" degree, tmp = "<<tmp<<", AngDist. = "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<", ThetaCmMin["<<reaction<<"] = "<<ThetaCmMin[reaction]/degree<<" degree"<<endl;
		  }
	      }
	    else
	      {
		data->Theta_cm = RandFlat::shoot(ThetaCmMin[reaction]/degree,ThetaCmMax[reaction]/degree)*degree;
		if(data->Theta_cm < ThetaCmMin[reaction] || ThetaCmMax[reaction] < data->Theta_cm)
		  {
		    cout<<"Error in GeneratePrimaryAction, reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" not possible (theta_cm = RandFlat::shoot("<<ThetaCmMin[reaction]<<","<<180.<<")*degree)"<<endl;
		    exit(15);
		  }
		tmp = RandFlat::shoot(0.,MaxAngularDistribution[reaction]);
		if(data->VerboseLevel > 4)
		  {
		    cout<<"GeneratePrimaryAction: reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" degree (RandFlat::shoot("<<ThetaCmMin[reaction]/degree<<","<<180.<<")), tmp = "<<tmp<<" (0-"<<MaxAngularDistribution[reaction]<<"), f(theta_cm) = "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<endl;
		  }
	      }
	    if(data->TestFlag)
	      {
		file<<data->Theta_cm/degree<<" "<<tmp<<" "<<AngularDistribution[reaction]->Eval(data->Theta_cm/degree)<<endl;
	      }
	  } while (tmp > AngularDistribution[reaction]->Eval(data->Theta_cm/degree));
	}//if(data->Transfer == 1)
      else
	{
	  //user provided angular distributions for elastic scattering
	  if(data->VerboseLevel > 4)
	    {
	      cout<<"shooting theta_cm from "<<ThetaCmMin[reaction]/radian<<" radian to "<<M_PI<<" radian (eex = "<<eex<<", reaction = "<<reaction<<", el.)"<<endl;
	    }
	  
	  do {
	    //old version from 0 to 180 degree
	    //theta_cm = ((log(data->OpticalFactor/data->OpticalExponent) - log(RandFlat::shoot(data->OpticalFactor/data->OpticalExponent,data->OpticalFactor/data->OpticalExponent*exp(-180.*data->OpticalExponent))))/data->OpticalExponent)*degree;
	    //new version from theta_cm_min to 180 degree, according to NR
	    //kwimmer 12.3.10 from thcmmin to thcmmax
	    data->Theta_cm = (-log(RandFlat::shoot(exp(-data->OpticalExponent*ThetaCmMin[reaction]/degree),exp(-data->OpticalExponent*ThetaCmMax[reaction]/degree)))/data->OpticalExponent)*degree;
	    if(data->Theta_cm < ThetaCmMin[reaction] || ThetaCmMax[reaction] < data->Theta_cm)
	      {
		cout<<"Error in GeneratePrimaryAction, reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" not possible (theta_cm = ((log("<<data->OpticalFactor<<"/"<<data->OpticalExponent<<") - log(RandFlat::shoot("<<data->OpticalFactor/data->OpticalExponent<<","<<data->OpticalFactor/data->OpticalExponent*exp(-180.*data->OpticalExponent)<<")))/"<<data->OpticalExponent<<")*degree"<<endl;
		exit(15);
	      }
	    tmp = RandFlat::shoot(0.,data->OpticalFactor*exp(-data->OpticalExponent*data->Theta_cm/degree));
	    if(data->VerboseLevel > 4)
	      {
		cout<<"GeneratePrimaryAction: reaction = "<<reaction<<", theta_cm = "<<data->Theta_cm/degree<<" degree (RandFlat::shoot("<<exp(-data->OpticalExponent*ThetaCmMin[reaction]/degree)<<","<<exp(-180.*data->OpticalExponent)<<")), tmp = "<<tmp<<" (0-"<<data->OpticalFactor*exp(-data->OpticalExponent*data->Theta_cm/degree)<<"), f(theta_cm) = "<<OpticalPotential[reaction - data->NumberOfLevel]->Eval(data->Theta_cm/degree)<<endl;
		cout<<"theta_cm = "<<data->Theta_cm/degree<<" degree, tmp = "<<tmp<<", AngDist. = "<<OpticalPotential[reaction - data->NumberOfLevel]->Eval(data->Theta_cm/degree)<<", ThetaCmMin["<<reaction<<"] = "<<ThetaCmMin[reaction]/degree<<" degree"<<endl;
	      }
	    if(data->TestFlag)
	      {
		file<<data->Theta_cm/degree<<" "<<tmp<<" "<<OpticalPotential[reaction - data->NumberOfLevel]->Eval(data->Theta_cm/degree)<<endl;
	      }
	  } while (tmp > OpticalPotential[reaction - data->NumberOfLevel]->Eval(data->Theta_cm/degree));
	}
    }//end of theta_cm
  //cout << data->Theta_cm/degree << endl;
  return;
}

void PrimaryGeneratorAction::ElasticTarget(int reaction, Element*& rc)
{
  //check which element from target material was chosen as "target"
  int i = reaction - data->NumberOfLevel;
  
  
  if(data->VerboseLevel > 4)
    cout<<"chose element "<<i<<" from "<<data->TargetMaterial->NumberOfElements()<<" in target ("<<data->TargetMaterial->Name()<<")"<<endl;
  
  rc = data->TargetMaterial->GetElement(i);
  
  if(rc == NULL)
    {
      cout<<"ERROR, tried to choose element "<<i<<" from "<<data->TargetMaterial->NumberOfElements()<<" in target ("<<data->TargetMaterial->Name()<<")"<<endl;
      exit(1);
    }
  
  return;
}


G4ParticleDefinition* PrimaryGeneratorAction::Definition(int Z, int N, double eex)
{
  if(Z+N > 4) //create ion from ion table
    return particleTable->GetIon(Z, Z+N, eex);
  else
    {
      if(Z == 1 && N == 0) //proton
	{
	  return G4Proton::ProtonDefinition();
	}
      else if(Z == 0 && N == 1) //proton
	{
	  return G4Neutron::NeutronDefinition();
	}
      else if(Z == 1 && N == 1) //deuteron
	{
	  return G4Deuteron::DeuteronDefinition();
	}
      else if(Z == 1 && N == 2) //triton
	{
	  return G4Triton::TritonDefinition();
	}
      else if(Z == 2 && N == 1) //3He
	{
	  return G4He3::He3Definition();
	}
      else if(Z == 2 && N == 2) //alpha
	{
	  return G4Alpha::AlphaDefinition();
	}
    }

  cerr<<"Error in "<<__PRETTY_FUNCTION__<<"shouldn't be able to reach this stage (Z = "<<Z<<", N = "<<N<<")"<<endl;
  exit(1);
}
